<p><b>xylar@lanl.gov</b> 2012-01-20 16:54:35 -0700 (Fri, 20 Jan 2012)</p><p>BRANCH COMMIT<br>
<br>
Main code:<br>
* Changed several variable names to be more descriptive (in consultation with Steve and Matt)<br>
unorm -> normalVelocity<br>
lsrf -> lowerSurface<br>
usrf -> upperSurface<br>
thck -> thickness<br>
thck_layer -> layerThickness<br>
temp -> temperature<br>
<br>
* Added nVertLevelsPlus2 = nVertLevels+2 as a dimension<br>
* gave temperature this vertical dimension<br>
<br>
* Added sigma coordinate fractions: layerThicknessFractions<br>
<br>
* Added test case 1: a dome with radius ~40 km and height ~700 m<br>
* Temperature is currently not being initialized in this test, but probably should be<br>
* Same with normalVelocity, layerThickness, bedTopography, lowerSurface and upperSurface<br>
<br>
* Added computation of layerThickness to the *very temporary* time step code<br>
<br>
Initial condition:<br>
* Added layerThicknessFractions and nVertLevelsPlus2 to icesheet/periodic_hex code<br>
* Added periodic_hex namelist for the periodic_hex dome test grid<br>
<br>
* Added a directory with a namelist for dome test runs under test_cases<br>
<br>
<br>
The land_ice core has now been tested, and works in that is reads in an input file,<br>
initializes a test case, performs a single (essentially empty) time step, and writes<br>
out both an output file and a restart file. The header of the restart file looks<br>
reasonable, but I haven't carefully verified the contents<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice/grid_gen/periodic_hex/Makefile
===================================================================
--- branches/land_ice/grid_gen/periodic_hex/Makefile 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/grid_gen/periodic_hex/Makefile 2012-01-20 23:54:35 UTC (rev 1406)
@@ -6,11 +6,11 @@
#LDFLAGS = -g -C
# pgf90
-FC = pgf90
-CC = pgcc
-FFLAGS = -r8 -O3
-CFLAGS = -O3
-LDFLAGS = -O3
+#FC = pgf90
+#CC = pgcc
+#FFLAGS = -r8 -O3
+#CFLAGS = -O3
+#LDFLAGS = -O3
# ifort
#FC = ifort
@@ -20,11 +20,11 @@
#LDFLAGS = -O3
# gfortran
-#FC = gfortran
-#CC = gcc
-#FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian
-#CFLAGS = -O3 -m64
-#LDFLAGS = -O3 -m64
+FC = gfortran
+CC = gcc
+FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian
+CFLAGS = -O3 -m64
+LDFLAGS = -O3 -m64
# absoft
#FC = f90
@@ -38,7 +38,7 @@
CPPFLAGS =
CPPINCLUDES =
INCLUDES = -I$(NETCDF)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf
+LIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdff
RM = rm -f
Modified: branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F
===================================================================
--- branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/icesheet/periodic_hex/module_write_netcdf.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -47,6 +47,7 @@
integer :: wrVarIDedgesOnVertex
integer :: wrVarIDcellsOnVertex
integer :: wrVarIDkiteAreasOnVertex
+ integer :: wrVarIDlayerThicknessFractions
integer :: wrVarIDthickness
integer :: wrVarIDbedTopography
integer :: wrVarIDbeta
@@ -211,6 +212,8 @@
dimlist( 1) = wrDimIDvertexDegree
dimlist( 2) = wrDimIDnVertices
nferr = nf_def_var(wr_ncid, 'kiteAreasOnVertex', NF_DOUBLE, 2, dimlist, wrVarIDkiteAreasOnVertex)
+ dimlist( 1) = wrDimIDnVertLevels
+ nferr = nf_def_var(wr_ncid, 'layerThicknessFractions', NF_DOUBLE, 1, dimlist, wrVarIDlayerThicknessFractions)
dimlist( 1) = wrDimIDnCells
dimlist( 2) = wrDimIDTime
nferr = nf_def_var(wr_ncid, 'thickness', NF_DOUBLE, 2, dimlist, wrVarIDthickness)
@@ -272,6 +275,7 @@
edgesOnVertex, &
cellsOnVertex, &
kiteAreasOnVertex, &
+ layerThicknessFractions, &
thickness, &
bedTopography, &
beta, &
@@ -319,6 +323,7 @@
integer, dimension(:,:), intent(in) :: edgesOnVertex
integer, dimension(:,:), intent(in) :: cellsOnVertex
real (kind=8), dimension(:,:), intent(in) :: kiteAreasOnVertex
+ real (kind=8), dimension(:), intent(in) :: layerThicknessFractions
real (kind=8), dimension(:,:), intent(in) :: thickness
real (kind=8), dimension(:,:), intent(in) :: bedTopography
real (kind=8), dimension(:,:), intent(in) :: beta
@@ -493,8 +498,12 @@
start2(2) = 1
count2( 1) = 3
count2( 2) = wrLocalnVertices
- nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start2, count2, kiteAreasOnVertex)
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDlayerThicknessFractions, start2, count2, kiteAreasOnVertex)
+ start1(1) = 1
+ count1( 1) = wrLocalnVertLevels
+ nferr = nf_put_vara_double(wr_ncid, wrVarIDkiteAreasOnVertex, start1, count1, layerThicknessFractions)
+
start2(2) = time
count2( 2) = wrLocalnCells
count2( 2) = 1
Modified: branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km
===================================================================
--- branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/icesheet/periodic_hex/namelist.input.domeTest.2km 2012-01-20 23:54:35 UTC (rev 1406)
@@ -1,6 +1,6 @@
&periodic_grid
nx = 30,
- ny = 30,
+ ny = 35,
dc = 2000.,
nVertLevels = 9,
nTracers = 1,
Modified: branches/land_ice/icesheet/periodic_hex/periodic_grid.F
===================================================================
--- branches/land_ice/icesheet/periodic_hex/periodic_grid.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/icesheet/periodic_hex/periodic_grid.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -23,6 +23,7 @@
real (kind=8), allocatable, dimension(:) :: latEdge, lonEdge, xEdge, yEdge, zEdge
real (kind=8), allocatable, dimension(:) :: latVertex, lonVertex, xVertex, yVertex, zVertex
real (kind=8), allocatable, dimension(:,:) :: weightsOnEdge, kiteAreasOnVertex
+ real (kind=8), allocatable, dimension(:) :: layerThicknessFractions
real (kind=8), allocatable, dimension(:,:) :: thickness, bedTopography, beta
real (kind=8), allocatable, dimension(:,:,:) :: normalVelocity
real (kind=8), allocatable, dimension(:,:,:,:) :: tracers
@@ -83,6 +84,8 @@
allocate(yVertex(nVertices))
allocate(zVertex(nVertices))
+ allocate(layerThicknessFractions(nVertLevels))
+
allocate(thickness(nCells,1))
allocate(bedTopography(nCells,1))
allocate(beta(nCells,1))
@@ -265,6 +268,9 @@
end do
+ ! evenly spaced layers for the sigma coordinate
+ layerThicknessFractions(:) = 1.0/real(nVertLevels)
+
!
! fill in initial conditions below
! NOTE: these initial conditions will likely be removed
@@ -304,6 +310,7 @@
edgesOnVertex, &
cellsOnVertex, &
kiteAreasOnVertex, &
+ layerThicknessFractions, &
thickness, &
bedTopography, &
beta, &
Modified: branches/land_ice/mpas/Makefile
===================================================================
--- branches/land_ice/mpas/Makefile 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/Makefile 2012-01-20 23:54:35 UTC (rev 1406)
@@ -241,7 +241,7 @@
CPPINCLUDES = -I../inc -I$(NETCDF)/include -I$(PAPI)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include -I$(PAPI)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf $(PAPILIBS)
+LIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdff $(PAPILIBS)
RM = rm -f
CPP = cpp -C -P -traditional
Modified: branches/land_ice/mpas/src/core_land_ice/Registry
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/Registry 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/Registry 2012-01-20 23:54:35 UTC (rev 1406)
@@ -12,17 +12,10 @@
namelist character land_ice_model config_run_duration none
namelist integer land_ice_model config_stats_interval 100
namelist logical land_ice_model config_h_ScaleWithMesh false
-%%%namelist real land_ice_model config_h_mom_eddy_visc2 0.0
-%%%namelist real land_ice_model config_h_mom_eddy_visc4 0.0
-%%%namelist real land_ice_model config_h_tracer_eddy_diff2 0.0
-%%%namelist real land_ice_model config_h_tracer_eddy_diff4 0.0
namelist integer land_ice_model config_thickness_adv_order 2
namelist integer land_ice_model config_tracer_adv_order 2
namelist logical land_ice_model config_positive_definite false
namelist logical land_ice_model config_monotonic false
-%%%namelist logical land_ice_model config_wind_stress false
-%%%namelist logical land_ice_model config_bottom_drag false
-%%%namelist real land_ice_model config_apvm_upwinding 0.5
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -46,6 +39,7 @@
dim TWENTYONE 21
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
+dim nVertLevelsPlus2 nVertLevels+2
dim nTracers nTracers
%
@@ -101,22 +95,11 @@
var persistent integer edgesOnVertex ( vertexDegree nVertices ) 0 iro edgesOnVertex mesh - -
var persistent integer cellsOnVertex ( vertexDegree nVertices ) 0 iro cellsOnVertex mesh - -
var persistent real kiteAreasOnVertex ( vertexDegree nVertices ) 0 iro kiteAreasOnVertex mesh - -
-%%%var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
-%%%var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
-%%%var persistent real fCell ( nCells ) 0 iro fCell mesh - -
-%%%var persistent real h_s ( nCells ) 0 iro h_s mesh - -
% Space needed for advection
var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
-% %% NOTE: the following arrays are needed to allow the use
-% %% of the module_advection.F w/o alteration
-% Space needed for deformation calculation weights
-var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
-var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
-var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
-
% Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
@@ -124,53 +107,35 @@
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
-%%%var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
+% Sigma coordinate: read from input, saved in restart and written to output
+var persistent real layerThicknessFractions ( nVertLevels ) 0 iro layerThicknessFractions mesh - -
+
% Prognostic variables: read from input, saved in restart, and written to output
-%%%var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
-%%%var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
-%%%var persistent real tracers ( nTracers nVertLevels nCells Time ) 2 iro tracers state - -
+var persistent real thickness ( nCells Time ) 2 iro thickness state - -
+var persistent real bedTopography ( nCells Time ) 2 iro bedTopography state - -
+% Note: beta should be moved to diagnostic fields when a process model comes online
+var persistent real beta ( nCells Time ) 2 iro beta state - -
+var persistent real temperature ( nVertLevelsPlus2 nCells Time ) 2 iro temperature state tracers dynamics
+var persistent real normalVelocity ( nVertLevels nEdges Time ) 2 iro normalVelocity state - -
% Tendency variables
-%%%var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
-var persistent real tend_thck_layer ( nVertLevels nCells Time ) 1 - thck_layer tend - -
-%%%var persistent real tend_tracers ( nTracers nVertLevels nCells Time ) 1 - tracers tend - -
+var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 - layerThickness tend - -
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
% Diagnostic fields: only written to output
-%%%var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
-%%%var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-%%%var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
-%%%var persistent real vorticity_cell ( nVertLevels nCells Time ) 2 o vorticity_cell state - -
-%%%var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
-var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
-%%%var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-%%%var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
-%%%var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
+var persistent real layerThickness ( nVertLevels nCells Time ) 2 o layerThickness state - -
+var persistent real lowerSurface ( nCells Time ) 2 o lowerSurface state - -
+var persistent real upperSurface ( nCells Time ) 2 o upperSurface state - -
% Other diagnostic variables: neither read nor written to any files
-%%%var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
-%%%var persistent real circulation ( nVertLevels nVertices Time ) 2 - circulation state - -
-%%%var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
-%%%var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
+var persistent real layerThicknessVertex ( nVertLevels nVertices Time ) 2 - layerThicknessVertex state - -
+var persistent real layerThicknessEdge ( nVertLevels nEdges Time ) 2 - layerThicknessEdge state - -
-% Prognostic variables for land-ice model: read from input, saved in restart, and written to output
-var persistent real thck ( nCells Time ) 2 iro thck state - -
-var persistent real lsrf ( nCells Time ) 2 iro lsrf state - -
-var persistent real beta ( nCells Time ) 2 iro beta state - -
-%%%var persistent real tracers ( nTracers nVertLevels nCells Time ) 2 iro tracers state - -
-var persistent real thck_layer ( nVertLevels nCells Time ) 2 iro thck_layer state - -
-var persistent real temp ( nVertLevels nCells Time ) 2 iro temp state tracers dynamics
-% whl - temp should be nVertLevels + 1?
-% Diagnostic fields for land-ice model: only written to output
-var persistent real unorm ( nVertLevels nEdges Time ) 2 iro unorm state - -
-var persistent real usrf ( nCells Time ) 2 iro usrf state - -
-
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -49,10 +49,10 @@
! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle
!!! real (kind=RKIND), dimension(:), pointer :: h_s, fCell, fEdge
-!!! real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, &
- real (kind=RKIND), dimension(:,:), pointer :: thck_layer, unorm, h_edge, &
+!!! real (kind=RKIND), dimension(:,:), pointer :: h, u, v, layerThicknessEdge, &
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness, normalVelocity, layerThicknessEdge, &
!!! pv_edge, pv_vertex, pv_cell, &
- h_vertex, weightsOnEdge
+ layerThicknessVertex, weightsOnEdge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -105,12 +105,12 @@
areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
weightsOnEdge => grid % weightsOnEdge % array
- thck_layer => state % thck_layer % array
- unorm => state % unorm % array
+ layerThickness => state % layerThickness % array
+ normalVelocity => state % normalVelocity % array
!!! v => state % v % array
tracers => state % tracers % array
- h_edge => state % h_edge % array
- h_vertex => state % h_vertex % array
+ layerThicknessEdge => state % layerThicknessEdge % array
+ layerThicknessVertex => state % layerThicknessVertex % array
!!! pv_edge => state % pv_edge % array
!!! pv_vertex => state % pv_vertex % array
!!! pv_cell => state % pv_cell % array
@@ -161,40 +161,40 @@
! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
do iLevel = 1,nVertLevels
! eg. GlobalFluidThickness top (Sum( h dA)) = Sum(cellVolume)
- cellVolume(iLevel,:) = thck_layer(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
+ cellVolume(iLevel,:) = layerThickness(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
cellArea(iLevel,:) = areaCell(1:nCellsSolve)
!!! volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
-!!! *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+!!! *layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
!!! volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
-!!! *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+!!! *pv_vertex(iLevel,1:nVerticesSolve)*layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ vertexVolume(iLevel,:) = layerThicknessVertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
!!! volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
-!!! *h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
+!!! *layerThicknessEdge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
!!! volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
!!! volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
-!!! refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(thck_layer(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
+!!! refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(layerThickness(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
do iEdge = 1,nEdgesSolve
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
!!! eoe = edgesOnEdge(j,iEdge)
!!! workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
-!!! q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
+!!! q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * layerThicknessEdge(iLevel,eoe)
end do
-!!! keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
+!!! keTend_CoriolisForce(iLevel,iEdge) = layerThicknessEdge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
iCell1 = cellsOnEdge(1,iEdge)
iCell2 = cellsOnEdge(2,iEdge)
-!!! refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
+!!! refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(layerThicknessEdge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
-!!! keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*h_edge(iLevel,iEdge)*u(iLevel,iEdge) &
+!!! keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*layerThicknessEdge(iLevel,iEdge)*u(iLevel,iEdge) &
!!! *gravity*(h(iLevel,iCell2)+h_s(iCell2) - h(iLevel,iCell1)-h_s(iCell1))/dcEdge(iEdge)
!!! peTend_DivThickness(iLevel,iCell1) = peTend_DivThickness(iLevel,iCell1) &
-!!! + h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+!!! + layerThicknessEdge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
!!! peTend_DivThickness(iLevel,iCell2) = peTend_DivThickness(iLevel,iCell2) &
-!!! - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+!!! - layerThicknessEdge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
end do
!!! peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_lifev.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -306,27 +306,27 @@
!
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:), pointer :: &
- thck, lsrf, beta
+ thickness, lowerSurface, beta
real (kind=RKIND), dimension(:,:), pointer :: &
- unorm
+ normalVelocity
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers
- integer :: index_temp
+ integer :: index_temperature
err = 0
- thck => state % thck % array
- lsrf => state % lsrf % array
+ thickness => state % thickness % array
+ lowerSurface => state % lowerSurface % array
beta => state % beta % array
tracers => state % tracers % array
- unorm => state % unorm % array
+ normalVelocity => state % normalVelocity % array
- index_temp = state % index_temp
+ index_temperature = state % index_temperature
- ! call lifev_velocity_solver_L1L2(thck, lsrf, beta, tracers(index_temp,:,:), unorm)
+ ! call lifev_velocity_solver_L1L2(thickness, lowerSurface, beta, tracers(index_temperature,:,:), normalVelocity)
!--------------------------------------------------------------------
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -136,7 +136,7 @@
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
- call mpas_reconstruct(mesh, block % state % time_levs(1) % state % unorm % array, &
+ call mpas_reconstruct(mesh, block % state % time_levs(1) % state % normalVelocity % array, &
block % state % time_levs(1) % state % uReconstructX % array, &
block % state % time_levs(1) % state % uReconstructY % array, &
block % state % time_levs(1) % state % uReconstructZ % array, &
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -41,6 +41,21 @@
!!! end do
!whl - removed remaining shallow water test cases
+ else if (config_test_case == 1) then
+ ! the dome test case
+ write(0,*) 'Setting up land ice test case 1'
+ write(0,*) ' -- Evolution of an isothermal dome'
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call land_ice_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
+ do i=2,nTimeLevs
+ call mpas_copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
else
write(0,*) 'Only test case 0 is currently supported.'
stop
@@ -49,6 +64,49 @@
end subroutine setup_land_ice_test_case
+ subroutine land_ice_test_case_1(grid, state)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup land ice test case 1: Isothermal dome
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+
+ real (kind=RKIND), parameter :: r0 = 60000.0*sqrt(0.125) !meters
+ real (kind=RKIND), parameter :: h0 = 2000.0*sqrt(0.125) !meters
+ real (kind=RKIND), parameter :: x0 = 30000.0 !meters
+ real (kind=RKIND), parameter :: y0 = 30000.0 !meters
+
+ integer :: iCell
+ real (kind=RKIND) :: r
+
+ real (kind=RKIND), dimension(:), pointer :: x, y, thickness
+
+ if(grid % on_a_sphere) then
+ print *, "Land ice test case 1 only supports planar grids."
+ stop
+ end if
+
+
+ x => grid % xCell % array
+ y => grid % yCell % array
+ thickness => state % thickness % array
+
+
+ do iCell = 1,grid % nCells
+ r = sqrt((x(iCell)-x0)**2 + (y(iCell)-y0)**2)
+ if(r < r0) then
+ thickness(iCell) = h0*sqrt(1 - (r/r0)**2)
+ else
+ thickness(iCell) = 0.0
+ end if
+ end do
+
+
+ end subroutine land_ice_test_case_1
+
!!! subroutine sw_test_case_1(grid, state)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-20 22:37:29 UTC (rev 1405)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-20 23:54:35 UTC (rev 1406)
@@ -27,8 +27,11 @@
type (block_type), pointer :: block
- integer :: err
+ real (kind=RKIND), dimension(:), pointer :: thickness, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ integer :: nVertLevels, nCells, iCell, iLevel, err
+
! if (trim(config_time_integration) == 'RK4') then
! call land_ice_rk4(domain, dt)
! else
@@ -40,6 +43,21 @@
! a temporary test that calls the velocity solver on the state at the current time
block => domain % blocklist
+ nCells = block % mesh % nCells
+ nVertLevels = block % mesh % nVertLevels
+ layerThicknessFractions => block % mesh % layerThicknessFractions % array
+
+ thickness => block % state % time_levs(1) % state % thickness % array
+ layerThickness => block % state % time_levs(1) % state % layerThickness % array
+
+ ! set up layerThickness from thickness using the layerThicknessFractions
+ do iCell = 1, nCells
+ do iLevel = 1, nVertLevels
+ layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
+ end do
+ end do
+
+ block => domain % blocklist
call land_ice_vel_solve(block % mesh, block % state % time_levs(1) % state, err)
block => domain % blocklist
@@ -90,15 +108,15 @@
block => domain % blocklist
do while (associated(block))
- block % state % time_levs(2) % state % unorm % array(:,:) = block % state % time_levs(1) % state % unorm % array(:,:)
+ block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
- block % state % time_levs(2) % state % thck_layer % array(:,:) = block % state % time_levs(1) % state % thck_layer % array(:,:)
+ block % state % time_levs(2) % state % layerThickness % array(:,:) = block % state % time_levs(1) % state % layerThickness % array(:,:)
do iCell=1,block % mesh % nCells ! couple tracers to h
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % thck_layer % array(k,iCell)
+ * block % state % time_levs(1) % state % layerThickness % array(k,iCell)
end do
end do
@@ -157,10 +175,10 @@
block => domain % blocklist
do while (associated(block))
-! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % unorm % array(:,:), &
+! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % normalVelocity % array(:,:), &
! block % mesh % nVertLevels, block % mesh % nEdges, &
! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % thck_layer % array(:,:), &
+ call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % layerThickness % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
@@ -174,21 +192,21 @@
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
-! provis % unorm % array(:,:) = block % state % time_levs(1) % state % unorm % array(:,:) &
-! + rk_substep_weights(rk_step) * block % tend % unorm % array(:,:)
- provis % thck_layer % array(:,:) = block % state % time_levs(1) % state % thck_layer % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % thck_layer % array(:,:)
+! provis % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:) &
+! + rk_substep_weights(rk_step) * block % tend % normalVelocity % array(:,:)
+ provis % layerThickness % array(:,:) = block % state % time_levs(1) % state % layerThickness % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % layerThickness % array(:,:)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
provis % tracers % array(:,k,iCell) = ( &
- block % state % time_levs(1) % state % thck_layer % array(k,iCell) * &
+ block % state % time_levs(1) % state % layerThickness % array(k,iCell) * &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % thck_layer % array(k,iCell)
+ ) / provis % layerThickness % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
- provis % unorm % array(:,:) = block % state % time_levs(1) % state % unorm % array(:,:)
+ provis % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
end if
call land_ice_compute_solve_diagnostics(dt, provis, block % mesh)
block => block % next
@@ -199,10 +217,10 @@
block => domain % blocklist
do while (associated(block))
-! block % state % time_levs(2) % state % unorm % array(:,:) = block % state % time_levs(2) % state % unorm % array(:,:) &
-! + rk_weights(rk_step) * block % tend % unorm % array(:,:)
- block % state % time_levs(2) % state % thck_layer % array(:,:) = block % state % time_levs(2) % state % thck_layer % array(:,:) &
- + rk_weights(rk_step) * block % tend % thck_layer % array(:,:)
+! block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(2) % state % normalVelocity % array(:,:) &
+! + rk_weights(rk_step) * block % tend % normalVelocity % array(:,:)
+ block % state % time_levs(2) % state % layerThickness % array(:,:) = block % state % time_levs(2) % state % layerThickness % array(:,:) &
+ + rk_weights(rk_step) * block % tend % layerThickness % array(:,:)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -228,17 +246,17 @@
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % state % time_levs(2) % state % thck_layer % array(k,iCell)
+ / block % state % time_levs(2) % state % layerThickness % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % state % time_levs(2) % state % unorm % array(:,:) = block % state % time_levs(1) % state % unorm % array(:,:)
+ block % state % time_levs(2) % state % normalVelocity % array(:,:) = block % state % time_levs(1) % state % normalVelocity % array(:,:)
end if
call land_ice_compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
- call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % unorm% array, &
+ call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % normalVelocity% array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
block % state % time_levs(2) % state % uReconstructZ % array, &
@@ -279,10 +297,10 @@
!!! real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge
!!! real (kind=RKIND), dimension(:,:), pointer :: vh
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, &
- h_edge, thck_layer, unorm, tend_thck_layer, &
-!!! h_edge, h, u, v, tend_h, tend_u, &
+ layerThicknessEdge, layerThickness, normalVelocity, tend_layerThickness, &
+!!! layerThicknessEdge, h, u, v, tend_h, tend_u, &
!!! circulation, vorticity, ke, pv_edge, divergence, &
- h_vertex
+ layerThicknessVertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r
@@ -300,10 +318,10 @@
!!! h => s % h % array
!!! u => s % u % array
!!! v => s % v % array
- thck_layer => s % thck_layer % array
- unorm => s % unorm % array
+ layerThickness => s % layerThickness % array
+ normalVelocity => s % normalVelocity % array
- h_edge => s % h_edge % array
+ layerThicknessEdge => s % layerThicknessEdge % array
!!! circulation => s % circulation % array
!!! vorticity => s % vorticity % array
!!! divergence => s % divergence % array
@@ -329,7 +347,7 @@
!!! fVertex => grid % fVertex % array
!!! fEdge => grid % fEdge % array
- tend_thck_layer => tend % thck_layer % array
+ tend_layerThickness => tend % layerThickness % array
! tend_u => tend % u % array
nCells = grid % nCells
@@ -346,19 +364,19 @@
!
! Compute height tendency for each cell
!
- tend_thck_layer(:,:) = 0.0
+ tend_layerThickness(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,nVertLevels
- flux = unorm(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_thck_layer(k,cell1) = tend_thck_layer(k,cell1) - flux
- tend_thck_layer(k,cell2) = tend_thck_layer(k,cell2) + flux
+ flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
+ tend_layerThickness(k,cell1) = tend_layerThickness(k,cell1) - flux
+ tend_layerThickness(k,cell2) = tend_layerThickness(k,cell2) + flux
end do
end do
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
- tend_thck_layer(k,iCell) = tend_thck_layer(k,iCell) / areaCell(iCell)
+ tend_layerThickness(k,iCell) = tend_layerThickness(k,iCell) / areaCell(iCell)
end do
end do
@@ -378,10 +396,10 @@
!!! do j = 1,nEdgesOnEdge(iEdge)
!!! eoe = edgesOnEdge(j,iEdge)
!!! workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
-!!! q = q + weightsOnEdge(j,iEdge) * unorm(k,eoe) * workpv * h_edge(k,eoe)
+!!! q = q + weightsOnEdge(j,iEdge) * normalVelocity(k,eoe) * workpv * layerThicknessEdge(k,eoe)
!!! end do
-!!! tend_unorm(k,iEdge) = &
+!!! tend_normalVelocity(k,iEdge) = &
!!! q &
!!! - ( ke(k,cell2) - ke(k,cell1) + &
!!! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
@@ -409,7 +427,7 @@
!!! workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
-!!! tend_unorm(k,iEdge) = workpv * vh(k,iEdge) - &
+!!! tend_normalVelocity(k,iEdge) = workpv * vh(k,iEdge) - &
!!! (ke(k,cell2) - ke(k,cell1) + &
!!! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
!!! ) / &
@@ -431,7 +449,7 @@
!!! u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
!!! -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
!!! u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
-!!! tend_unorm(k,iEdge) = tend_unorm(k,iEdge) + u_diffusion
+!!! tend_normalVelocity(k,iEdge) = tend_normalVelocity(k,iEdge) + u_diffusion
!!! end do
!!!! end do
!!! end if
@@ -459,7 +477,7 @@
!!! do k=1,nVertLevels
-!!! delsq_unorm(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+!!! delsq_normalVelocity(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
!!! -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
!!! end do
@@ -535,7 +553,7 @@
!!! if(config_wind_stress) then
!!! do iEdge=1,grid % nEdges
!!! tend_u(1,iEdge) = tend_u(1,iEdge) &
-!!! + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+!!! + u_src(1,iEdge)/rho_ref/layerThicknessEdge(1,iEdge)
!!! end do
!!! endif
@@ -548,8 +566,8 @@
!!! + ke(1,cellsOnEdge(2,iEdge)))
!!! tend_u(1,iEdge) = tend_u(1,iEdge) &
-!!! - 1.0e-3*unorm(1,iEdge) &
-!!! *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+!!! - 1.0e-3*normalVelocity(1,iEdge) &
+!!! *sqrt(2.0*ke_edge)/layerThicknessEdge(1,iEdge)
!!! end do
!!! endif
@@ -583,10 +601,10 @@
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND), dimension(:,:), pointer :: unorm, h_edge
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThicknessEdge
- unorm => s % unorm % array
- h_edge => s % h_edge % array
+ normalVelocity => s % normalVelocity % array
+ layerThicknessEdge => s % layerThicknessEdge % array
dcEdge => grid % dcEdge % array
deriv_two => grid % deriv_two % array
dvEdge => grid % dvEdge % array
@@ -613,7 +631,7 @@
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = unorm(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * tracer_edge
tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
end do
@@ -658,14 +676,14 @@
endif
!-- if u > 0:
- if (unorm(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * unorm(k,iEdge) * h_edge(k,iEdge) * ( &
+ if (normalVelocity(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
!-- else u <= 0:
else
- flux = dvEdge(iEdge) * unorm(k,iEdge) * h_edge(k,iEdge) * ( &
+ flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
@@ -715,7 +733,7 @@
endif
- flux = dvEdge(iEdge) * unorm(k,iEdge) * h_edge(k,iEdge) * ( &
+ flux = dvEdge(iEdge) * normalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge) * ( &
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
@@ -754,7 +772,7 @@
!!! *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
-!!! flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
+!!! flux = dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
!!! tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
!!! tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
!!! end do
@@ -791,9 +809,9 @@
!!! do k=1,grid % nVertLevels
!!! do iTracer=1, grid % nTracers
!!! delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
-!!! + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+!!! + dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
!!! delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
-!!! - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+!!! - dvEdge(iEdge) * layerThicknessEdge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
!!! end do
!!! end do
@@ -857,12 +875,12 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, areaTriangle
!!! real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge
!!! real (kind=RKIND), dimension(:,:), pointer :: vh
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, h_edge, &
- thck_layer, unorm, tend_thck_layer, &
-!!! real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, layerThicknessEdge, &
+ layerThickness, normalVelocity, tend_layerThickness, &
+!!! real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, layerThicknessEdge, h, u, v, tend_h, tend_u, &
!!! circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, &
!!! gradPVn, gradPVt, divergence, vorticity_cell, &
- h_vertex
+ layerThicknessVertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, &
edgesOnCell, edgesOnEdge, edgesOnVertex, &
boundaryEdge, boundaryCell
@@ -872,13 +890,13 @@
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
- thck_layer => s % thck_layer % array
- unorm => s % unorm % array
+ layerThickness => s % layerThickness % array
+ normalVelocity => s % normalVelocity % array
!!! v => s % v % array
!!! vh => s % vh % array
- h_edge => s % h_edge % array
- h_vertex => s % h_vertex % array
- tend_thck_layer => s % thck_layer % array
+ layerThicknessEdge => s % layerThicknessEdge % array
+ layerThicknessVertex => s % layerThicknessVertex % array
+ tend_layerThickness => s % layerThickness % array
!!! tend_u => s % u % array
!!! circulation => s % circulation % array
!!! vorticity => s % vorticity % array
@@ -935,7 +953,7 @@
!
! Compute height on cell edges at velocity locations
- ! Namelist options control the order of accuracy of the reconstructed h_edge value
+ ! Namelist options control the order of accuracy of the reconstructed layerThicknessEdge value
!
coef_3rd_order = 0.
@@ -949,7 +967,7 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
- h_edge(k,iEdge) = 0.5 * (thck_layer(k,cell1) + thck_layer(k,cell2))
+ layerThicknessEdge(k,iEdge) = 0.5 * (layerThickness(k,cell1) + layerThickness(k,cell2))
end do
end if
end do
@@ -971,33 +989,33 @@
!-- if not a boundary cell
if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * thck_layer(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * thck_layer(k,cell2)
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * layerThickness(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * layerThickness(k,cell2)
!-- all edges of cell 1
do i=1, grid % nEdgesOnCell % array (cell1)
d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * thck_layer(k,grid % CellsOnCell % array (i,cell1))
+ deriv_two(i+1,1,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell1))
end do
!-- all edges of cell 2
do i=1, grid % nEdgesOnCell % array (cell2)
d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * thck_layer(k,grid % CellsOnCell % array (i,cell2))
+ deriv_two(i+1,2,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell2))
end do
endif
!-- if u > 0:
- if (unorm(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(thck_layer(k,cell1) + thck_layer(k,cell2)) &
+ if (normalVelocity(k,iEdge) > 0) then
+ layerThicknessEdge(k,iEdge) = &
+ 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else unorm <= 0:
+ !-- else normalVelocity <= 0:
else
- h_edge(k,iEdge) = &
- 0.5*(thck_layer(k,cell1) + thck_layer(k,cell2)) &
+ layerThicknessEdge(k,iEdge) = &
+ 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
end if
@@ -1023,25 +1041,25 @@
!-- if not a boundary cell
if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * thck_layer(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * thck_layer(k,cell2)
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * layerThickness(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * layerThickness(k,cell2)
!-- all edges of cell 1
do i=1, grid % nEdgesOnCell % array (cell1)
d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * thck_layer(k,grid % CellsOnCell % array (i,cell1))
+ deriv_two(i+1,1,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell1))
end do
!-- all edges of cell 2
do i=1, grid % nEdgesOnCell % array (cell2)
d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * thck_layer(k,grid % CellsOnCell % array (i,cell2))
+ deriv_two(i+1,2,iEdge) * layerThickness(k,grid % CellsOnCell % array (i,cell2))
end do
endif
- h_edge(k,iEdge) = &
- 0.5*(thck_layer(k,cell1) + thck_layer(k,cell2)) &
+ layerThicknessEdge(k,iEdge) = &
+ 0.5*(layerThickness(k,cell1) + layerThickness(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
end do ! do k
@@ -1054,7 +1072,7 @@
! 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
!
- unorm(:,nEdges+1) = 0.0
+ normalVelocity(:,nEdges+1) = 0.0
!
! Compute circulation and relative vorticity at each vertex
@@ -1062,8 +1080,8 @@
!!! circulation(:,:) = 0.0
!!! do iEdge=1,nEdges
!!! do k=1,nVertLevels
-!!! circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * unorm(k,iEdge)
-!!! circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * unorm(k,iEdge)
+!!! circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * normalVelocity(k,iEdge)
+!!! circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * normalVelocity(k,iEdge)
!!! end do
!!! end do
!!! do iVertex=1,nVertices
@@ -1082,12 +1100,12 @@
!!! cell2 = cellsOnEdge(2,iEdge)
!!! if (cell1 <= nCells) then
!!! do k=1,nVertLevels
-!!! divergence(k,cell1) = divergence(k,cell1) + unorm(k,iEdge)*dvEdge(iEdge)
+!!! divergence(k,cell1) = divergence(k,cell1) + normalVelocity(k,iEdge)*dvEdge(iEdge)
!!! enddo
!!! endif
!!! if(cell2 <= nCells) then
!!! do k=1,nVertLevels
-!!! divergence(k,cell2) = divergence(k,cell2) - unorm(k,iEdge)*dvEdge(iEdge)
+!!! divergence(k,cell2) = divergence(k,cell2) - normalVelocity(k,iEdge)*dvEdge(iEdge)
!!! enddo
!!! end if
!!! end do
@@ -1106,7 +1124,7 @@
!!! do i=1,nEdgesOnCell(iCell)
!!! iEdge = edgesOnCell(i,iCell)
!!! do k=1,nVertLevels
-!!! ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * unorm(k,iEdge)**2.0
+!!! ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * normalVelocity(k,iEdge)**2.0
!!! end do
!!! end do
!!! do k=1,nVertLevels
@@ -1122,7 +1140,7 @@
!!! do i=1,nEdgesOnEdge(iEdge)
!!! eoe = edgesOnEdge(i,iEdge)
!!! do k = 1,nVertLevels
-!!! v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * unorm(k, eoe)
+!!! v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * normalVelocity(k, eoe)
!!! end do
!!! end do
!!! end do
@@ -1136,7 +1154,7 @@
!!! do j=1,nEdgesOnEdge(iEdge)
!!! eoe = edgesOnEdge(j,iEdge)
!!! do k=1,nVertLevels
-!!! vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * unorm(k,eoe) * h_edge(k,eoe)
+!!! vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * normalVelocity(k,eoe) * layerThicknessEdge(k,eoe)
!!! end do
!!! end do
!!! end do
@@ -1149,13 +1167,13 @@
!
do iVertex = 1,nVertices
do k=1,nVertLevels
- h_vertex(k,iVertex) = 0.0
+ layerThicknessVertex(k,iVertex) = 0.0
do i=1,grid % vertexDegree
- h_vertex(k,iVertex) = h_vertex(k,iVertex) + thck_layer(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+ layerThicknessVertex(k,iVertex) = layerThicknessVertex(k,iVertex) + layerThickness(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
- h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
+ layerThicknessVertex(k,iVertex) = layerThicknessVertex(k,iVertex) / areaTriangle(iVertex)
-!!! pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
+!!! pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / layerThicknessVertex(k,iVertex)
end do
end do
@@ -1233,12 +1251,12 @@
!
!!! do iEdge = 1,nEdges
!!! do k = 1,nVertLevels
-!!! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * unorm(k,iEdge) * dt * gradPVn(k,iEdge)
+!!! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * normalVelocity(k,iEdge) * dt * gradPVn(k,iEdge)
!!! enddo
!!! enddo
!
- ! set pv_edge = fEdge / h_edge at boundary points
+ ! set pv_edge = fEdge / layerThicknessEdge at boundary points
!
! if (maxval(boundaryEdge).ge.0) then
! do iEdge = 1,nEdges
@@ -1248,13 +1266,13 @@
! if(boundaryEdge(k,iEdge).eq.1) then
! v(k,iEdge) = 0.0
! if(cell1.gt.0) then
- ! h1 = thck_layer(k,cell1)
+ ! h1 = layerThickness(k,cell1)
! pv_edge(k,iEdge) = fEdge(iEdge) / h1
- ! h_edge(k,iEdge) = h1
+ ! layerThicknessEdge(k,iEdge) = h1
! else
- ! h2 = thck_layer(k,cell2)
+ ! h2 = layerThickness(k,cell2)
! pv_edge(k,iEdge) = fEdge(iEdge) / h2
- ! h_edge(k,iEdge) = h2
+ ! layerThicknessEdge(k,iEdge) = h2
! endif
! endif
! enddo
Added: branches/land_ice/test_cases/land_ice/dome/namelist.input
===================================================================
--- branches/land_ice/test_cases/land_ice/dome/namelist.input (rev 0)
+++ branches/land_ice/test_cases/land_ice/dome/namelist.input 2012-01-20 23:54:35 UTC (rev 1406)
@@ -0,0 +1,27 @@
+&land_ice_model
+ config_test_case = 1
+ config_time_integration = 'RK4'
+ config_dt = 60.0
+ config_start_time = '0000-01-01_00:00:00'
+ config_stop_time = '0000-01-01_00:01:00'
+ config_run_duration = '00:01:00'
+ config_stats_interval = 0
+ config_h_ScaleWithMesh = .false.
+ config_thickness_adv_order = 2
+ config_tracer_adv_order = 2
+ config_positive_definite = .false.
+ config_monotonic = .false.
+/
+
+&io
+ config_input_name = 'grid.nc'
+ config_output_name = 'output.nc'
+ config_restart_name = 'restart.nc'
+ config_output_interval = '00_00:01:00'
+ config_frames_per_outfile = 0
+/
+
+&restart
+ config_restart_interval = '00_00:01:00'
+ config_do_restart = .false.
+/
</font>
</pre>