<p><b>mpetersen@lanl.gov</b> 2010-05-07 13:17:39 -0600 (Fri, 07 May 2010)</p><p>Merged trunk changes, including the lateral boundary condition, to the z-level branch. Tested with isopycnal 3-layer double gyre domain.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -2,7 +2,7 @@
MODEL_FORMULATION = -DLANL_FORMULATION
#EXPAND_LEVELS = -DEXPAND_LEVELS=26
-#FILE_OFFSET = -DOFFSET64BIT
+FILE_OFFSET = -DOFFSET64BIT
#########################
# Section for Zoltan TPL
@@ -75,14 +75,37 @@
"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
Modified: branches/ocean_projects/z_level_mrp/mpas/namelist.input.hyd_atmos
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/namelist.input.hyd_atmos 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/namelist.input.hyd_atmos 2010-05-07 19:17:39 UTC (rev 260)
@@ -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/ocean_projects/z_level_mrp/mpas/namelist.input.ocean
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/namelist.input.ocean 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/namelist.input.ocean 2010-05-07 19:17:39 UTC (rev 260)
@@ -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/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -19,10 +19,9 @@
mpas_interface.o: module_advection.o module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry 2010-05-07 19:17:39 UTC (rev 260)
@@ -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
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -96,7 +96,7 @@
do_the_cell = .true.
do i=1,n
- if (cell_list(i) <= 0) do_the_cell = .false.
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
end do
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_test_cases.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_test_cases.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -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/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -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)
@@ -671,10 +671,10 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_divergence(nVertLevels, nCells))
- allocate(delsq_u(nVertLevels, nEdges))
- allocate(delsq_circulation(nVertLevels, nVertices))
- allocate(delsq_vorticity(nVertLevels, nVertices))
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
@@ -849,7 +849,7 @@
if ( h_theta_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -1519,9 +1519,9 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( num_scalars, grid % nEdges) :: h_flux
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( num_scalars, grid % nEdges+1) :: h_flux
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: scale_out, scale_in
real (kind=RKIND), dimension( num_scalars ) :: s_max, s_min, s_max_update, s_min_update
integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -19,10 +19,9 @@
mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-07 19:17:39 UTC (rev 260)
@@ -92,11 +92,12 @@
var integer maxLevelsCell ( nCells ) iro kmaxCell - -
var integer maxLevelsEdge ( nEdges ) iro kmaxEdge - -
var real hZLevel ( nVertLevels ) iro hZLevel - -
-var real zmidZLevel ( nVertLevels ) iro zmidZLevel - -
-var real zbotZLevel ( nVertLevels ) iro zbotZLevel - -
+var real zMidZLevel ( nVertLevels ) iro zMidZLevel - -
+var real zBotZLevel ( nVertLevels ) iro zBotZLevel - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
var real u_src ( nVertLevels nEdges ) iro u_src - -
# Prognostic variables: read from input, saved in restart, and written to output
@@ -120,21 +121,18 @@
var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
-# mrp 100112:
-var real zmid ( nVertLevels nCells Time ) o zmid - -
-var real zbot ( nVertLevels nCells Time ) o zbot - -
+var real zMid ( nVertLevels nCells Time ) o zMid - -
+var real zBot ( nVertLevels nCells Time ) o zBot - -
var real zSurface ( nCells Time ) o zSurface - -
-var real zmidEdge ( nVertLevels nEdges Time ) o zmidEdge - -
-var real zbotEdge ( nVertLevels nEdges Time ) o zbotEdge - -
+var real zMidEdge ( nVertLevels nEdges Time ) o zMidEdge - -
+var real zBotEdge ( nVertLevels nEdges Time ) o zBotEdge - -
var real zSurfaceEdge ( nEdges Time ) o zSurfaceEdge - -
-var real pmid ( nVertLevels nCells Time ) o pmid - -
-var real pbot ( nVertLevels nCells Time ) o pbot - -
-var real pmidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
+var real pMid ( nVertLevels nCells Time ) o pMid - -
+var real pBot ( nVertLevels nCells Time ) o pBot - -
+var real pMidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
var real w ( nVertLevels nCells Time ) o w - -
-# mrp 100112 end
-
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
var real circulation ( nVertLevels nVertices Time ) - circulation - -
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -21,14 +21,18 @@
type (domain_type), intent(inout) :: domain
- integer :: i, iCell, iEdge, iVtx, iLevel, iTracer
+ integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
+
+ ! mrp 100507: for diagnostic output
+ integer :: iTracer
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
hZLevel, zmidZLevel, zbotZLevel
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
integer :: nCells, nEdges, nVertices, nVertLevels
+ ! mrp 100507 end: for diagnostic output
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -79,164 +83,16 @@
stop
end if
- ! mrp 100121:
- !
- ! Initialize u_src, rho, alpha
- ! This is a temporary fix until everything is specified correctly
- ! in an initial conditions file.
- !
block_ptr => domain % blocklist
do while (associated(block_ptr))
- h => block_ptr % time_levs(1) % state % h % array
- u => block_ptr % time_levs(1) % state % u % array
- rho => block_ptr % time_levs(1) % state % rho % array
- tracers => block_ptr % time_levs(1) % state % tracers % array
- u_src => block_ptr % mesh % u_src % array
- xCell => block_ptr % mesh % xCell % array
- yCell => block_ptr % mesh % yCell % array
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % time_levs(1) % state, &
+ block_ptr % time_levs(i) % state)
+ end do
- hZLevel => block_ptr % mesh % hZLevel % array
- zmidZLevel => block_ptr % mesh % zmidZLevel % array
- zbotZLevel => block_ptr % mesh % zbotZLevel % array
-
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
- ! Momentum forcing u_src
- if (config_test_case > 0) then
- ! for shallow water test cases:
- u_src = 0.0
- elseif (config_test_case == 0 ) then
- ! for rectangular basin:
- !do iEdge=1,nEdges
- ! u_src(1:nVertLevels, iEdge) = u_src(1, iEdge) / nVertLevels
- !end do
- endif
-
- ! set rho directly:
- !delta_rho=0.0
- !do iLevel = 1,nVertLevels
- ! rho(iLevel,1) = delta_rho*(iLevel-1)
- !enddo
- !rho(:,1) = rho(:,1) - sum(rho(:,1))/nVertLevels + 1000.0
- !do iLevel = 1,nVertLevels
- ! rho(iLevel,:) = rho(iLevel,1)
- !enddo
-
- tracers = 0.0
- do iCell = 1,nCells
- ! for three layer test
-! tracers(index_salinity,1,iCell) = 2.0 ! salinity
-! tracers(index_salinity,2,iCell) = 6.0 ! salinity
-! tracers(index_salinity,3,iCell) = 13.0 ! salinity
- do iLevel = 1,nVertLevels
-! if (xCell(iCell)<500*1e3.and.ycell(iCell)<1000*1e3) then
-! tracers(index_tracer1,iLevel,iCell) = 1.0
-! tracers(index_tracer2,iLevel,iCell) = 1.0e4
-! endif
-! if (ycell(iCell)>1400*1e3.and.ycell(iCell)<1500*1e3) then
-! tracers(index_tracer1,iLevel,iCell) = 1.0
-! tracers(index_tracer2,iLevel,iCell) = 1.0e4
-! endif
-! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
-! tracers(index_temperature,iLevel,iCell) = 6.0-iLevel ! temperature
-! tracers(index_temperature,iLevel,iCell) = 6.0-iLevel+yCell(iCell)/4000.e3 ! temperature
-! tracers(index_salinity,iLevel,iCell) = 35.0 ! salinity
-
- ! for 20 layer test
- tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
-
-! tracers(index_temperature,iLevel,iCell) = &
-! 10.0* xCell(iCell)/2500.e3
-! tracers(index_salinity,iLevel,iCell) = &
-! 34.0 + 2* yCell(iCell)/4000.e3
- tracers(index_tracer1,iLevel,iCell) = 1.0
- tracers(index_tracer2,iLevel,iCell) = &
- (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
-! tracers(index_tracer1,iLevel,iCell) = mod(ilevel,2)
- rho(iLevel,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(index_temperature,iLevel,iCell) &
- + 7.6e-4*tracers(index_salinity,iLevel,iCell))
- enddo
- enddo
-
-#ifdef EXPAND_LEVELS
- print '(10a)', 'EXPAND_LEVELS compiler flag is on.', &
- ' Copying h and u from k=1 to other levels.'
- print '(a,i5)', 'EXPAND_LEVELS =', EXPAND_LEVELS
- print '(a,i5)', 'nVertLevels =', nVertLevels
- do iCell=1,nCells
- ! make the total thickness equal to the single-layer thickness:
- h(1:nVertLevels, iCell) = h(1,iCell) / nVertLevels
- end do
-
- do iEdge=1,nEdges
- u(2:nVertLevels, iEdge) = u(1,iEdge)
- end do
-#else
- print '(10a)', 'EXPAND_LEVELS compiler flag is off.'
-#endif
-
-
- ! mrp 100426: initialize z-level variables.
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- hZLevel = h(:,1)
- zmidZLevel(1) = -0.5*hZLevel(1)
- zbotZLevel(1) = -hZLevel(1)
- do iLevel = 2,nVertLevels
- zmidZLevel(iLevel) = zbotZLevel(iLevel-1)-0.5*hZLevel(iLevel)
- zbotZLevel(iLevel) = zbotZLevel(iLevel-1)- hZLevel(iLevel)
- enddo
- if (config_vert_grid_type.eq.'isopycnal') then
- print *, ' Using isopycnal coordinates'
- elseif (config_vert_grid_type.eq.'zlevel') then
- print *, ' Using z-level coordinates'
- else
- print *, ' Incorrect choice of config_vert_grid_type:',&
- config_vert_grid_type
- stop
- endif
-
- ! mrp 100426 end: initialize z-level variables.
-
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, &
- block_ptr % time_levs(i) % state)
- end do
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min u_src max u_src ', &
- ' min h max h ',&
- ' hZLevel zmidZlevel', &
- ' zbotZlevel'
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,:)), maxval(u(iLevel,:)), &
- minval(u_src(iLevel,:)), maxval(u_src(iLevel,:)), &
- minval(h(iLevel,:)), maxval(h(iLevel,:)), &
- hZLevel(iLevel),zmidZlevel(iLevel),zbotZlevel(iLevel)
- enddo
-
- ! print some diagnostics
- print '(10a)', 'itracer ilevel min tracer max tracer'
- do iTracer=1,num_tracers
- do iLevel = 1,nVertLevels
- print '(2i5,20es12.4)', iTracer,ilevel, &
- minval(tracers(itracer,iLevel,:)), maxval(tracers(itracer,iLevel,:))
- enddo
- enddo
-
- block_ptr => block_ptr % next
+ block_ptr => block_ptr % next
end do
- ! mrp 100121 end
end subroutine setup_sw_test_case
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -118,7 +118,6 @@
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
-
! --- update halos for diagnostic variables
block => domain % blocklist
@@ -133,10 +132,9 @@
block => domain % blocklist
do while (associated(block))
-
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -210,6 +208,7 @@
+ rk_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell)
end do
end do
+
block => block % next
end do
@@ -281,8 +280,8 @@
type (grid_meta), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wHat, rho0Inv
real (kind=RKIND), dimension(:), pointer :: &
@@ -291,7 +290,8 @@
real (kind=RKIND), dimension(:,:), pointer :: &
vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, &
tend_h, tend_u, Fv, circulation, vorticity, ke, pv_edge, &
- divergence, MontPot, pmidZLevel, zMidEdge, zbotEdge
+ divergence, MontPot, pmidZLevel, zMidEdge, zBotEdge
+
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
@@ -299,10 +299,8 @@
real (kind=RKIND) :: u_diffusion, visc, dist
real (kind=RKIND), dimension(:), allocatable:: fluxVert
-!ocean
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
-!ocean
visc = config_visc
@@ -317,13 +315,11 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
vh => s % vh % array
- !mrp 100112:
MontPot => s % MontPot % array
pmidZLevel => s % pmidZLevel % array
- zmidEdge => s % zmidEdge % array
- zbotEdge => s % zbotEdge % array
+ zBotEdge => s % zBotEdge % array
zSurfaceEdge=> s % zSurfaceEdge % array
- !mrp 100112 end
+ zMidEdge => s % zMidEdge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -351,9 +347,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
-!ocean
u_src => grid % u_src % array
-!ocean
!
@@ -364,13 +358,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
end do
end if
- if (cell2 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -446,7 +440,6 @@
!
tend_u(:,:) = 0.0
rho0Inv = 1.0/config_rho0
- allocate(fluxVert(0:nVertLevels))
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -510,7 +503,10 @@
- ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
end do
+ end do
+ do iEdge=1,grid % nEdgesSolve
+
! forcing in top layer only
tend_u(1,iEdge) = tend_u(1,iEdge) + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
@@ -519,7 +515,11 @@
tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- u(nVertLevels,iEdge)/(100.0*86400.0)
- ! vertical mixing
+ enddo
+
+ ! vertical mixing
+ allocate(fluxVert(0:nVertLevels))
+ do iEdge=1,grid % nEdgesSolve
fluxVert(0) = 0.0
fluxVert(nVertLevels) = 0.0
do k=nVertLevels-1,1,-1
@@ -538,9 +538,9 @@
tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
enddo
-
end do
deallocate(fluxVert)
+
#endif
#ifdef NCAR_FORMULATION
@@ -591,7 +591,6 @@
nEdges, nCells, nVertLevels
real (kind=RKIND) :: flux, tracer_edge, r, &
Tmid(num_tracers,grid % nVertLevels)
-
real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
@@ -630,7 +629,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
do iTracer=1,num_tracers
tracer_edge = 0.5 * ( tracers(iTracer,k,cell1) &
@@ -649,7 +648,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
if (u(k,iEdge)>0.0) then
upwindCell = cell1
@@ -668,7 +667,6 @@
endif
-
! mrp 100409 computation of h tendency was, only had div.(huT) term
!do iCell=1,grid % nCellsSolve
! do k=1,grid % nVertLevelsSolve
@@ -791,10 +789,11 @@
! leave it separate for now for clarity.
! default kappa is 2e3m^2/s, from 2e7cm^2/s in the POP Users Guide p. 39
allocate(tr_flux(num_tracers,nVertLevels,nEdges))
+ tr_flux(:,:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
do iTracer=1,num_tracers
@@ -812,7 +811,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
do iTracer=1,num_tracers
tr_div(iTracer,k,cell1) = tr_div(iTracer,k,cell1) &
@@ -820,7 +819,7 @@
enddo
enddo
endif
- if(cell2 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
do iTracer=1,num_tracers
tr_div(iTracer,k,cell2) = tr_div(iTracer,k,cell2) &
@@ -843,13 +842,15 @@
deallocate(tr_flux, tr_div)
! mrp 100501 end: Horizontal tracer diffusion
- ! print some diagnostics
- !do iTracer=1,num_tracers
- !do k = 1,nVertLevels
- ! print '(2i5,20es12.4)', iTracer,k, &
- ! minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))
- !enddo
- !enddo
+ ! print some diagnostics - for debugging
+! print *, 'after vertical mixing',&
+! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
+! do iTracer=1,num_tracers
+! do k = 1,nVertLevels
+! print '(2i5,20es12.4)', iTracer,k, &
+! minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))
+! enddo
+! enddo
end subroutine compute_scalar_tend
@@ -875,6 +876,8 @@
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, pbot_temp
integer :: nCells, nEdges, nVertices, nVertLevels
+
+
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zSurface, hZLevel, zSurfaceEdge
@@ -886,7 +889,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -908,7 +911,6 @@
pv_cell => s % pv_cell % array
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
- !mrp 100112:
rho => s % rho % array
tracers => s % tracers % array
zmid => s % zmid % array
@@ -921,7 +923,6 @@
MontPot => s % MontPot % array
zSurface => s % zSurface % array
zSurfaceEdge=> s % zSurfaceEdge % array
- !mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -947,7 +948,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -955,24 +956,39 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
+ elseif(cell1 <= nCells) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = h(k,cell1)
+ end do
+ elseif(cell2 <= nCells) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = h(k,cell2)
+ end do
end if
end do
+
!
+ ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
+ ! used to when reading for edges that do not exist
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
end do
end if
- if (verticesOnEdge(2,iEdge) > 0) then
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -984,6 +1000,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -991,12 +1008,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -1044,7 +1061,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -1068,14 +1085,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -1083,14 +1099,12 @@
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
-
+
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -1102,7 +1116,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -1111,16 +1124,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -1131,7 +1142,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -1140,31 +1150,29 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
+
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
@@ -1174,29 +1182,11 @@
enddo
!
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- h1 = 0.0
- h2 = 0.0
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
-
- !mrp 100324:
- !
! equation of state
!
+ ! For a isopycnal model, density should remain constant.
+ if (config_vert_grid_type.eq.'zlevel') then
+
do iCell=1,nCells
do k=1,nVertLevels
! Linear equation of state, for the time being
@@ -1205,11 +1195,11 @@
+ 7.6e-4*tracers(index_salinity,k,iCell))
end do
end do
+
+ endif
!mrp 100324 end
- !mrp 100112:
- !
! For Isopycnal model.
! Compute mid- and bottom-depth of each layer, from bottom up.
! Then compute mid- and bottom-pressure of each layer, and
@@ -1220,14 +1210,14 @@
! h_s is ocean topography: elevation above lowest point,
! and is positive. z coordinates are pos upward, and z=0
! is at lowest ocean point.
- zbot(nVertLevels,iCell) = h_s(iCell)
- zmid(nVertLevels,iCell) = zbot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
+ zBot(nVertLevels,iCell) = h_s(iCell)
+ zMid(nVertLevels,iCell) = zBot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
do k=nVertLevels-1,1,-1
- zbot(k,iCell) = zbot(k+1,iCell) + h(k+1,iCell)
- zmid(k,iCell) = zbot(k,iCell) + 0.5*h(k,iCell)
+ zBot(k,iCell) = zBot(k+1,iCell) + h(k+1,iCell)
+ zMid(k,iCell) = zBot(k,iCell) + 0.5*h(k,iCell)
end do
- ! rather than using zbot(0,iCell), I am adding another 2D variable.
- zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
+ ! rather than using zBot(0,iCell), I am adding another 2D variable.
+ zSurface(iCell) = zBot(1,iCell) + h(1,iCell)
! assume atmospheric pressure at the surface is zero for now.
pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell)
@@ -1244,12 +1234,11 @@
end do
end do
- !mrp 100112 end
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if(cell1.gt.0.and.cell2.gt.0) then
+ if(cell1<=nCells .and. cell2<=nCells) then
do k=1,nVertLevels
zBotEdge(k,iEdge) = (zBot(k,cell1)+zBot(k,cell2))/2.0
zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
@@ -1285,13 +1274,13 @@
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1300,7 +1289,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -1310,21 +1299,21 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -16,10 +16,9 @@
mpas_interface.o: module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Registry 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_sw/Registry 2010-05-07 19:17:39 UTC (rev 260)
@@ -81,7 +81,8 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
# Prognostic variables: read from input, saved in restart, and written to output
var real u ( nVertLevels nEdges Time ) iro u - -
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_sw/module_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_sw/module_test_cases.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_sw/module_test_cases.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -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**2.0 * 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/ocean_projects/z_level_mrp/mpas/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_sw/module_time_integration.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_sw/module_time_integration.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -125,7 +125,7 @@
do while (associated(block))
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -297,13 +297,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
end do
end if
- if (cell2 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -343,6 +343,7 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
+
tend_u(k,iEdge) = &
q &
+ u_diffusion &
@@ -404,7 +405,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
@@ -450,7 +451,7 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -495,7 +496,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -503,7 +504,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
@@ -511,16 +512,22 @@
end do
!
+ ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
+ ! used to when reading for edges that do not exist
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
end do
end if
- if (verticesOnEdge(2,iEdge) > 0) then
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -532,6 +539,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -539,12 +547,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -580,7 +588,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -604,14 +612,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -623,10 +630,8 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -638,7 +643,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -647,16 +651,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -667,7 +669,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -676,30 +677,28 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
+
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
! Modify PV edge with upstream bias.
!
@@ -712,32 +711,38 @@
!
! set pv_edge = fEdge / h_edge at boundary points
!
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
+ ! if (maxval(boundaryEdge).ge.0) then
+ ! do iEdge = 1,nEdges
+ ! cell1 = cellsOnEdge(1,iEdge)
+ ! cell2 = cellsOnEdge(2,iEdge)
+ ! do k = 1,nVertLevels
+ ! if(boundaryEdge(k,iEdge).eq.1) then
+ ! v(k,iEdge) = 0.0
+ ! if(cell1.gt.0) then
+ ! h1 = h(k,cell1)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
+ ! h_edge(k,iEdge) = h1
+ ! else
+ ! h2 = h(k,cell2)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
+ ! h_edge(k,iEdge) = h2
+ ! endif
+ ! endif
+ ! enddo
+ ! enddo
+ ! endif
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -746,7 +751,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -756,22 +761,22 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
- tend_u => tend % u % array
+ boundaryEdge => grid % boundaryEdge % array
+ tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: branches/ocean_projects/z_level_mrp/mpas/src/driver/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/driver/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/driver/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -10,10 +10,9 @@
mpas.o: module_subdriver.o
clean:
- $(RM) *.o *.mod
+ $(RM) *.o *.mod *.f90
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../core_$(CORE)
- $(RM) $*.f90
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -35,13 +35,12 @@
module_io_output.o: module_grid_types.o module_dmpar.o module_sort.o module_configure.o
clean:
- $(RM) *.o *.mod libframework.a
+ $(RM) *.o *.mod *.f90 libframework.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES)
- $(RM) $*.f90
.c.o:
$(CC) $(CFLAGS) $(CPPFLAGS) $(CPPINCLUDES) -c $<
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_block_decomp.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_block_decomp.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_block_decomp.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -142,7 +142,11 @@
edgeIDListLocal(:) = edgeIDList(:)
do i=1,nEdges
- if (hash_search(h, cellsOnEdge(1,i))) then
+ do j=1,maxCells
+ if (cellsOnEdge(j,i) /= 0) exit
+ end do
+if (j > maxCells) write(0,*) 'Error in block_decomp_partitioned_edge_list: edge/vertex is not adjacent to any valid cells'
+ if (hash_search(h, cellsOnEdge(j,i))) then
lastEdge = lastEdge + 1
edgeIDList(lastEdge) = edgeIDListLocal(i)
else
@@ -231,8 +235,10 @@
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ end if
end if
end do
end do
@@ -264,10 +270,12 @@
write(0,*) 'block_decomp_add_halo: Somehow we don''t have the right number of non-ghost cells.'
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
- local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
- k = k + 1
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
+ k = k + 1
+ end if
end if
end do
end do
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_dmpar.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_dmpar.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_dmpar.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -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
@@ -688,8 +707,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- integer, dimension(nOwnedList), intent(in) :: arrayIn
- integer, dimension(nNeededList), intent(inout) :: arrayOut
+ integer, dimension(*), intent(in) :: arrayIn
+ integer, dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -765,7 +784,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -778,8 +797,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- integer, dimension(dim1,nOwnedList), intent(in) :: arrayIn
- integer, dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ integer, dimension(dim1,*), intent(in) :: arrayIn
+ integer, dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -857,7 +876,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -869,8 +888,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- real (kind=RKIND), dimension(nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -946,7 +965,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -959,8 +978,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1038,7 +1057,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -1051,8 +1070,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,dim2,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,dim2,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1130,7 +1149,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:,:) = arrayIn(:,:,:)
+ arrayOut(:,:,1:nNeededList) = arrayIn(:,:,1:nOwnedList)
end if
#endif
@@ -1142,7 +1161,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- integer, dimension(nField), intent(in) :: field
+ integer, dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1169,7 +1188,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- integer, dimension(ds:de,1:nField), intent(in) :: field
+ integer, dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1203,7 +1222,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(nField), intent(in) :: field
+ real (kind=RKIND), dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1230,7 +1249,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1264,7 +1283,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1302,7 +1321,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- integer, dimension(nField), intent(inout) :: field
+ integer, dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1329,7 +1348,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- integer, dimension(ds:de,1:nField), intent(inout) :: field
+ integer, dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1358,7 +1377,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(nField), intent(inout) :: field
+ real (kind=RKIND), dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1385,7 +1404,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1415,7 +1434,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1449,7 +1468,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1
- real (kind=RKIND), dimension(dim1), intent(inout) :: array
+ real (kind=RKIND), dimension(*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1509,7 +1528,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2
- real (kind=RKIND), dimension(dim1,dim2), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1573,7 +1592,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, dim3
- real (kind=RKIND), dimension(dim1,dim2,dim3), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -21,6 +21,7 @@
interface io_input_field
+ module procedure io_input_field0dReal
module procedure io_input_field1dReal
module procedure io_input_field2dReal
module procedure io_input_field3dReal
@@ -804,7 +805,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnCell % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -812,7 +814,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnCell % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -820,7 +823,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnCell % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
end if
end do
@@ -834,7 +838,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnEdge % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -842,7 +847,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnEdge % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
end if
end do
@@ -854,7 +860,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnEdge % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
end if
end do
@@ -868,7 +875,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnVertex % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -876,7 +884,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnVertex % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
end if
end do
@@ -993,10 +1002,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
@@ -1031,6 +1040,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"
@@ -1057,6 +1077,33 @@
end subroutine io_input_get_dimension
+ subroutine io_input_field0dReal(input_obj, field)
+
+ implicit none
+
+ type (io_input_object), intent(in) :: input_obj
+ type (field0dReal), intent(inout) :: field
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+ integer :: varID
+ integer, dimension(1) :: start1, count1
+
+ start1(1) = 1
+ count1(1) = 1
+
+#include "input_field0dreal.inc"
+
+#if (RKIND == 8)
+ nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#else
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#endif
+
+ end subroutine io_input_field0dReal
+
+
subroutine io_input_field1dReal(input_obj, field)
implicit none
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -27,6 +27,7 @@
interface io_output_field
+ module procedure io_output_field0dReal
module procedure io_output_field1dReal
module procedure io_output_field2dReal
module procedure io_output_field3dReal
@@ -126,9 +127,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
@@ -187,8 +188,12 @@
domain % blocklist % mesh % edgesOnEdge % array(j,i))
end do
do j=domain % blocklist % mesh % nEdgesOnEdge % array(i)+1,2*domain % blocklist % mesh % maxEdges
- edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
+ if(domain % blocklist % mesh % nEdgesOnEdge % array(i) .eq. 0) then
+ edgesOnEdge(j,i) = domain % blocklist % mesh % nEdgesSolve + 1
+ else
+ edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
domain % blocklist % mesh % nEdgesOnEdge % array(i))
+ endif
end do
end do
do i=1,domain % blocklist % mesh % nVerticesSolve
@@ -336,6 +341,35 @@
end subroutine io_output_init
+ subroutine io_output_field0dReal(output_obj, field)
+
+ implicit none
+
+ type (io_output_object), intent(in) :: output_obj
+ type (field0dReal), intent(inout) :: field
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+ integer :: varID
+ integer, dimension(1) :: start1, count1
+
+ start1(1) = 1
+ count1(1) = 1
+
+#include "output_field0dreal.inc"
+
+#if (RKIND == 8)
+ nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#else
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#endif
+
+ nferr = nf_sync(output_obj % wr_ncid)
+
+ end subroutine io_output_field0dReal
+
+
subroutine io_output_field1dReal(output_obj, field)
implicit none
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_zoltan_interface.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_zoltan_interface.F 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_zoltan_interface.F 2010-05-07 19:17:39 UTC (rev 260)
@@ -1,9 +1,10 @@
module zoltan_interface
use zoltan
- use mpi
implicit none
+ include 'mpif.h'
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Data for reordering cells
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: branches/ocean_projects/z_level_mrp/mpas/src/operators/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/operators/Makefile 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/operators/Makefile 2010-05-07 19:17:39 UTC (rev 260)
@@ -10,10 +10,9 @@
module_vector_reconstruction.o:
clean:
- $(RM) *.o *.mod libops.a
+ $(RM) *.o *.mod *.f90 libops.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
- $(RM) $*.f90
Modified: branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c 2010-05-07 14:00:22 UTC (rev 259)
+++ branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c 2010-05-07 19:17:39 UTC (rev 260)
@@ -387,10 +387,20 @@
fortprintf(fd, " allocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(g %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -414,16 +424,28 @@
else {
fortprintf(fd, " allocate(g %% %s)</font>
<font color="black">", var_ptr->name_in_code);
fortprintf(fd, " allocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " allocate(g %% %s %% array(", var_ptr->name_in_code);
- dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
- dimlist_ptr = dimlist_ptr->next;
- while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (var_ptr->ndims > 0) {
+ fortprintf(fd, " allocate(g %% %s %% array(", var_ptr->name_in_code);
+ dimlist_ptr = var_ptr->dimlist;
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
+ while (dimlist_ptr) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ dimlist_ptr = dimlist_ptr->next;
+ }
+ fortprintf(fd, "))</font>
<font color="red">");
+
}
- fortprintf(fd, "))</font>
<font color="red">");
-
if (var_ptr->iostreams & INPUT0)
fortprintf(fd, " g %% %s %% ioinfo %% input = .true.</font>
<font color="gray">", var_ptr->name_in_code);
else
@@ -473,9 +495,15 @@
fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="red">", var_ptr2->super_array);
}
else {
- fortprintf(fd, " deallocate(g %% %s %% array)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="blue">", var_ptr->name_in_code);
+ if (var_ptr->ndims > 0) {
+ fortprintf(fd, " deallocate(g %% %s %% array)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="blue">", var_ptr->name_in_code);
+ }
+ else {
+ fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="gray">", var_ptr->name_in_code);
+ }
var_ptr = var_ptr->next;
}
}
@@ -508,10 +536,20 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(s %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -538,10 +576,20 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
fortprintf(fd, " allocate(s %% %s %% array(", var_ptr->name_in_code);
dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -1026,7 +1074,10 @@
fortprintf(fd, " deallocate(super_%s%id)</font>
<font color="red">", vtype, var_ptr->ndims);
}
else {
- fortprintf(fd, " block %% time_levs(1) %% state %% %s %% scalar = %s%id %% scalar</font>
<font color="blue">", var_ptr->name_in_code, vtype, var_ptr->ndims);
+ if (var_ptr->timedim)
+ fortprintf(fd, " block %% time_levs(1) %% state %% %s %% scalar = %s%id %% scalar</font>
<font color="blue">", var_ptr->name_in_code, vtype, var_ptr->ndims);
+ else
+ fortprintf(fd, " block %% mesh %% %s %% scalar = %s%id %% scalar</font>
<font color="black">", var_ptr->name_in_code, vtype, var_ptr->ndims);
}
fortprintf(fd, " end if</font>
<font color="black"></font>
<font color="gray">");
@@ -1086,10 +1137,10 @@
/*
- * Generate code to read 1d, 2d, 3d time-invariant fields
+ * Generate code to read 0d, 1d, 2d, 3d time-invariant fields
*/
for(j=0; j<2; j++) {
- for(i=1; i<=3; i++) {
+ for(i=0; i<=3; i++) {
if (j == 0) {
sprintf(fname, "input_field%idinteger.inc", i);
ivtype = INTEGER;
@@ -1497,7 +1548,10 @@
}
else {
fortprintf(fd, " %s%id %% ioinfo %% fieldName = \'%s\'</font>
<font color="red">", vtype, var_ptr->ndims, var_ptr->name_in_file);
- fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% time_levs(1) %% state %% %s %% scalar</font>
<font color="blue">", vtype, var_ptr->ndims, var_ptr->name_in_code);
+ if (var_ptr->timedim)
+ fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% time_levs(1) %% state %% %s %% scalar</font>
<font color="blue">", vtype, var_ptr->ndims, var_ptr->name_in_code);
+ else
+ fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% mesh %% %s %% scalar</font>
<font color="gray">", vtype, var_ptr->ndims, var_ptr->name_in_code);
}
if (var_ptr->timedim)
@@ -1518,10 +1572,10 @@
/*
- * Generate code to write 1d, 2d, 3d time-invariant fields
+ * Generate code to write 0d, 1d, 2d, 3d time-invariant fields
*/
for(j=0; j<2; j++) {
- for(i=1; i<=3; i++) {
+ for(i=0; i<=3; i++) {
if (j == 0) {
sprintf(fname, "output_field%idinteger.inc", i);
ivtype = INTEGER;
</font>
</pre>