<p><b>xylar@lanl.gov</b> 2011-05-10 14:41:41 -0600 (Tue, 10 May 2011)</p><p>Merged the trunk to the land_ice branch<br>
Copied the Registry and module_time_integration.F from core_sw to core_land_ice<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/land_ice
___________________________________________________________________
Added: svn:mergeinfo
+ /trunk:618-825
Copied: branches/land_ice/documents/MPASProjectDesignTemplate.tex (from rev 825, trunk/documents/MPASProjectDesignTemplate.tex)
===================================================================
--- branches/land_ice/documents/MPASProjectDesignTemplate.tex (rev 0)
+++ branches/land_ice/documents/MPASProjectDesignTemplate.tex 2011-05-10 20:41:41 UTC (rev 826)
@@ -0,0 +1,82 @@
+\documentclass[11pt]{report}
+
+\usepackage{epsf,amsmath,amsfonts}
+\usepackage{graphicx}
+
+\begin{document}
+
+\title{Title: \\
+Requirements and Design}
+\author{MPAS Development Team}
+
+\maketitle
+\tableofcontents
+
+%-----------------------------------------------------------------------
+
+\chapter{Summary}
+
+The purpose of this section is to summarize what capability is to be added to the MPAS system through this design process. It should be clear what new code will do that the current code does not. Summarizing the primary challenges with respect to software design and implementation is also appropriate for this section. Finally, this statement should contain general statement with regard to what is ``success.''
+
+%figure template
+%\begin{figure}
+% \center{\includegraphics[width=14cm]{./Figure1.pdf}}
+% \caption{A graphical representation of the discrete boundary.}
+% \label{Figure1}
+%\end{figure}
+
+
+
+
+%-----------------------------------------------------------------------
+
+\chapter{Requirements}
+
+\section{Requirement: XXX}
+Date last modified: 2011/01/05 \\
+Contributors: (add your name to this list if it does not appear) \\
+
+Each requirement is to be listed under a ''section'' heading, as there will be a one-to-one correspondence between requirements, design, proposed implementation and testing.
+
+Requirements should not discuss technical software issues, but rather focus on model capability. To the extent possible, requirements should be relatively independent of each other, thus allowing a clean design solution, implementation and testing plan.
+
+%-----------------------------------------------------------------------
+
+\chapter{Algorithmic Formulations}
+
+\section{Design Solution: XXX}
+Date last modified: 2011/01/05 \\
+Contributors: (add your name to this list if it does not appear) \\
+
+For each requirement, there is a design solution that is intended to meet that requirement. Design solutions can include detailed technical discussions of PDEs, algorithms, solvers and similar, as well as technical discussion of performance issues. In general, this section should steer away from a detailed discussion of low-level software issues such as variable declarations, interfaces and sequencing.
+
+
+%-----------------------------------------------------------------------
+
+\chapter{Design and Implementation}
+
+\section{Implementation: XXX}
+Date last modified: 2011/01/05 \\
+Contributors: (add your name to this list if it does not appear) \\
+
+This section should detail the plan for implementing the design solution for requirement XXX. In general, this section is software-centric with a focus on software implementation. Pseudo code is appropriate in this section. Links to actual source code are appropriate. Project management items, such as svn branches, timelines and staffing are also appropriate.
+
+How do we typeset pseudo code?
+\begin{verbatim}
+ verbatim?
+\end{verbatim}
+
+%-----------------------------------------------------------------------
+
+\chapter{Testing}
+
+\section{Testing and Validation: XXX}
+Date last modified: 2011/01/05 \\
+Contributors: (add your name to this list if it does not appear) \\
+
+How will XXX be tested? i.e. how will be we know when we have met requirement XXX. Will these unit tests be included in the ongoing going forward?
+
+
+%-----------------------------------------------------------------------
+
+\end{document}
Property changes on: branches/land_ice/mpas
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/ocean_projects/vert_adv_mrp:704-745
/trunk/mpas:618-825
Modified: branches/land_ice/mpas/Makefile
===================================================================
--- branches/land_ice/mpas/Makefile 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/Makefile 2011-05-10 20:41:41 UTC (rev 826)
@@ -38,7 +38,7 @@
"CC = cc" \
"SFC = ftn" \
"SCC = cc" \
- "FFLAGS = -i4 -r8 -gopt -O2 -Mvect=nosse -Kieee" \
+ "FFLAGS = -i4 -r8 -gopt -O2 -Mvect=nosse -Kieee -convert big_endian" \
"CFLAGS = -fast" \
"LDFLAGS = " \
"CORE = $(CORE)" \
@@ -50,7 +50,7 @@
"CC = mpicc" \
"SFC = pgf90" \
"SCC = pgcc" \
- "FFLAGS = -r8 -O3" \
+ "FFLAGS = -r8 -O3 -byteswapio" \
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
@@ -62,7 +62,7 @@
"CC = pgcc" \
"SFC = pgf90" \
"SCC = pgcc" \
- "FFLAGS = -i4 -r8 -g -O2" \
+ "FFLAGS = -i4 -r8 -g -O2 -byteswapio" \
"CFLAGS = -fast" \
"LDFLAGS = " \
"CORE = $(CORE)" \
@@ -74,7 +74,7 @@
"CC = pgcc" \
"SFC = pgf90" \
"SCC = pgcc" \
- "FFLAGS = -r8 -O0 -g -Mbounds -Mchkptr" \
+ "FFLAGS = -r8 -O0 -g -Mbounds -Mchkptr -byteswapio" \
"CFLAGS = -O0 -g" \
"LDFLAGS = -O0 -g -Mbounds -Mchkptr" \
"CORE = $(CORE)" \
@@ -86,7 +86,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -O3" \
+ "FFLAGS = -real-size 64 -O3 -convert big_endian" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
@@ -98,7 +98,7 @@
"CC = mpicc" \
"SFC = gfortran" \
"SCC = gcc" \
- "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
@@ -110,7 +110,7 @@
"CC = gcc" \
"SFC = gfortran" \
"SCC = gcc" \
- "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8 -fconvert=big-endian" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3 -m64" \
"CORE = $(CORE)" \
@@ -122,7 +122,7 @@
"CC = mpicc" \
"SFC = g95" \
"SCC = gcc" \
- "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8 -fendian=big" \
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
@@ -134,7 +134,7 @@
"CC = gcc" \
"SFC = g95" \
"SCC = gcc" \
- "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+ "FFLAGS = -O3 -ffree-line-length-huge -r8 -fendian=big" \
"CFLAGS = -O3" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
Modified: branches/land_ice/mpas/namelist.input.ocean
===================================================================
--- branches/land_ice/mpas/namelist.input.ocean 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/namelist.input.ocean 2011-05-10 20:41:41 UTC (rev 826)
@@ -41,8 +41,12 @@
config_vert_viscosity = 1.0e-4
config_vert_diffusion = 1.0e-4
/
+&eos
+ config_eos_type = 'linear'
+/
&advection
- config_vert_tracer_adv = 'upwind'
+ config_vert_tracer_adv = 'spline'
+ config_vert_tracer_adv_order = 3
config_tracer_adv_order = 2
config_thickness_adv_order = 2
config_positive_definite = .false.
Modified: branches/land_ice/mpas/namelist.input.sw
===================================================================
--- branches/land_ice/mpas/namelist.input.sw 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/namelist.input.sw 2011-05-10 20:41:41 UTC (rev 826)
@@ -13,6 +13,8 @@
config_tracer_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
+ config_wind_stress = .false.
+ config_bottom_drag = .false.
/
&io
Modified: branches/land_ice/mpas/src/core_hyd_atmos/Registry
===================================================================
--- branches/land_ice/mpas/src/core_hyd_atmos/Registry 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_hyd_atmos/Registry 2011-05-10 20:41:41 UTC (rev 826)
@@ -23,6 +23,7 @@
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
+namelist character io config_decomp_file_prefix graph.info.part.
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
Modified: branches/land_ice/mpas/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_hyd_atmos/module_time_integration.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_hyd_atmos/module_time_integration.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -1937,10 +1937,7 @@
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells )
!
- VTX_LOOP: do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
- end do
+ do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex = 0.0
do i=1,grid % vertexDegree
@@ -1950,7 +1947,7 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
- end do VTX_LOOP
+ end do
! tdr
Modified: branches/land_ice/mpas/src/core_land_ice/Registry
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/Registry 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_land_ice/Registry 2011-05-10 20:41:41 UTC (rev 826)
@@ -15,10 +15,13 @@
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
+namelist logical sw_model config_wind_stress false
+namelist logical sw_model config_bottom_drag false
namelist real sw_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
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
+namelist character io config_decomp_file_prefix graph.info.part.
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -111,6 +114,7 @@
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 - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
Modified: branches/land_ice/mpas/src/core_land_ice/module_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/module_time_integration.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_land_ice/module_time_integration.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -269,8 +269,10 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+ real (kind=RKIND) :: ke_edge
-
h => s % h % array
u => s % u % array
v => s % v % array
@@ -308,26 +310,20 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ u_src => grid % u_src % array
!
! Compute height tendency for each cell
!
tend_h(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- 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 <= 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
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
@@ -504,6 +500,28 @@
end if
+ ! Compute u (velocity) tendency from wind stress (u_src)
+ 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)
+ end do
+ endif
+
+ if (config_bottom_drag) then
+ do iEdge=1,grid % nEdges
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+ ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
+ + ke(1,cellsOnEdge(2,iEdge)))
+
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ - 1.0e-3*u(1,iEdge) &
+ *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+ end do
+ endif
+
end subroutine compute_tend
@@ -1092,10 +1110,7 @@
! 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) > nCells) cycle VTX_LOOP
- end do
+ do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex(k,iVertex) = 0.0
do i=1,grid % vertexDegree
@@ -1105,7 +1120,7 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
end do
- end do VTX_LOOP
+ end do
!
Modified: branches/land_ice/mpas/src/core_ocean/Registry
===================================================================
--- branches/land_ice/mpas/src/core_ocean/Registry 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_ocean/Registry 2011-05-10 20:41:41 UTC (rev 826)
@@ -7,9 +7,10 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist integer sw_model config_stats_interval 100
-namelist character io config_input_name grid.nc
-namelist character io config_output_name output.nc
-namelist character io config_restart_name restart.nc
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
+namelist character io config_decomp_file_prefix graph.info.part.
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -30,7 +31,9 @@
namelist real vmix config_vmixTanhDiffMin 1.0e-5
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
-namelist character advection config_vert_tracer_adv 'centered'
+namelist character eos config_eos_type linear
+namelist character advection config_vert_tracer_adv stencil
+namelist integer advection config_vert_tracer_adv_order 4
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
namelist logical advection config_positive_definite false
@@ -106,7 +109,7 @@
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 real deriv_two ( FIFTEEN TWO nEdges ) 0 - deriv_two mesh - -
var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
# !! NOTE: the following arrays are needed to allow the use
@@ -120,11 +123,17 @@
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
# Arrays for z-level version of mpas-ocean
-var persistent integer maxLevelsCell ( nCells ) 0 iro kmaxCell mesh - -
-var persistent integer maxLevelsEdge ( nEdges ) 0 iro kmaxEdge mesh - -
+var persistent integer maxLevelCell ( nCells ) 0 iro maxLevelCell mesh - -
+var persistent integer maxLevelEdgeTop ( nEdges ) 0 - maxLevelEdgeTop mesh - -
+var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
+var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
+var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
-var persistent real zMidZLevel ( nVertLevels ) 0 iro zMidZLevel mesh - -
-var persistent real zTopZLevel ( nVertLevelsP1 ) 0 iro zTopZLevel mesh - -
+var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
+var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
+var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
+var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
+var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
@@ -150,31 +159,24 @@
var persistent real tend_tracer2 ( nVertLevels nCells Time ) 1 - tracer2 tend tracers testing
# Diagnostic fields: only written to output
-var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
+var persistent real v ( nVertLevels nEdges Time ) 2 - 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 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 h_vertex ( nVertLevels nVertices Time ) 2 o h_vertex state - -
+var persistent real vorticity ( nVertLevels nVertices Time ) 2 - vorticity state - -
+var persistent real pv_edge ( nVertLevels nEdges Time ) 2 - pv_edge state - -
+var persistent real h_edge ( nVertLevels nEdges Time ) 2 - h_edge state - -
+var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real ke_edge ( nVertLevels nEdges Time ) 2 o ke_edge state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
+var persistent real ke_edge ( nVertLevels nEdges Time ) 2 - ke_edge state - -
+var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 - 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 zMid ( nVertLevels nCells Time ) 2 o zMid state - -
-var persistent real zTop ( nVertLevelsP1 nCells Time ) 2 o zTop state - -
-var persistent real zMidEdge ( nVertLevels nEdges Time ) 2 o zMidEdge state - -
-var persistent real zTopEdge ( nVertLevelsP1 nEdges Time ) 2 o zTopEdge state - -
-var persistent real p ( nVertLevels nCells Time ) 2 o p state - -
-var persistent real pTop ( nVertLevelsP1 nCells Time ) 2 o pTop state - -
-var persistent real pZLevel ( nVertLevels nCells Time ) 2 o pZLevel state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
+var persistent real pressure ( nVertLevels nCells Time ) 2 o pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
-var persistent real ssh ( nCells Time ) 2 o ssh state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -182,7 +184,6 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-# xsad 10-02-05:
# Globally reduced diagnostic variables: only written to output
var persistent real areaCellGlobal ( Time ) 2 o areaCellGlobal state - -
var persistent real areaEdgeGlobal ( Time ) 2 o areaEdgeGlobal state - -
@@ -191,5 +192,4 @@
var persistent real volumeCellGlobal ( Time ) 2 o volumeCellGlobal state - -
var persistent real volumeEdgeGlobal ( Time ) 2 o volumeEdgeGlobal state - -
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
-# xsad 10-02-05 end
Modified: branches/land_ice/mpas/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/land_ice/mpas/src/core_ocean/module_global_diagnostics.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_ocean/module_global_diagnostics.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -37,7 +37,7 @@
real (kind=RKIND) :: areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, pv_edge, pv_vertex, &
- pv_cell, gradPVn, gradPVt, zMid, zTop, p, pTop, MontPot, wTop, rho, tracerTemp
+ pv_cell, gradPVn, gradPVt, pressure, MontPot, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
@@ -80,11 +80,8 @@
pv_cell => state % pv_cell % array
gradPVn => state % gradPVn % array
gradPVt => state % gradPVt % array
- zMid => state % zMid % array
- zTop => state % zTop % array
- p => state % p % array
- pTop => state % pTop % array
MontPot => state % MontPot % array
+ pressure => state % pressure % array
variableIndex = 0
! h
@@ -156,30 +153,12 @@
gradPVt(:,1:nEdgesSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! zMid
+ ! pressure
variableIndex = variableIndex + 1
- call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zMid(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
- ! zTop
- variableIndex = variableIndex + 1
- call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
- ! p
- variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- p(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ pressure(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! pTop
- variableIndex = variableIndex + 1
- call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- pTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
! MontPot
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
@@ -309,22 +288,10 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeEdgeGlobal
- ! zMid
+ ! pressure
variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
-
- ! zTop
- variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
-
- ! p
- variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
- ! pTop
- variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
-
! MontPot
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
Modified: branches/land_ice/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/land_ice/mpas/src/core_ocean/module_mpas_core.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_ocean/module_mpas_core.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -1,6 +1,8 @@
module mpas_core
use mpas_framework
+ use dmpar
+ use test_cases
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -13,7 +15,6 @@
use configure
use grid_types
- use test_cases
implicit none
@@ -21,10 +22,23 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
+ type (dm_info) :: dminfo
-
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ 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'
+ call init_ZLevel(domain)
+ else
+ print *, ' Incorrect choice of config_vert_grid_type:',&
+ config_vert_grid_type
+ call dmpar_abort(dminfo)
+ endif
+
+ call compute_maxLevel(domain)
+
!
! Initialize core
!
@@ -62,15 +76,45 @@
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
+ integer :: i, iEdge, iCell, k
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-
+
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ ! initialize velocities and tracers on land to be -1e34
+ ! The reconstructed velocity on land will have values not exactly
+ ! -1e34 due to the interpolation of reconstruction.
+
+ do iEdge=1,block % mesh % nEdges
+ ! mrp 101115 note: in order to include flux boundary conditions, the following
+ ! line will need to change. Right now, set boundary edges between land and
+ ! water to have zero velocity.
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeTop % array(iEdge)+1 &
+ :block % mesh % maxLevelEdgeBot % array(iEdge), iEdge) = 0.0
+
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeBot % array(iEdge)+1: &
+ block % mesh % nVertLevels,iEdge) = -1e34
+ end do
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(1) % state % tracers % array( &
+ :, block % mesh % maxLevelCell % array(iCell)+1 &
+ :block % mesh % nVertLevels,iCell) = -1e34
+ end do
- if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
+ if (.not. config_do_restart) then
+ block % state % time_levs(1) % state % xtime % scalar = 0.0
+ else
+ do i=2,nTimeLevs
+ call copy_state(block % state % time_levs(i) % state, &
+ block % state % time_levs(1) % state)
+ end do
+ endif
end subroutine mpas_init_block
@@ -108,14 +152,18 @@
! Move time level 2 fields back into time level 1 for next time step
call shift_time_levels_state(domain % blocklist % state)
- if (mod(itimestep, config_output_interval) == 0) then
- call write_output_frame(output_obj, output_frame, domain)
- end if
- if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
- if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
- call output_state_for_domain(restart_obj, domain, restart_frame)
- restart_frame = restart_frame + 1
- end if
+ if (config_output_interval > 0) then
+ if (mod(itimestep, config_output_interval) == 0) then
+ call write_output_frame(output_obj, output_frame, domain)
+ end if
+ endif
+ if (config_restart_interval > 0) then
+ if (mod(itimestep, config_restart_interval) == 0) then
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
+ call output_state_for_domain(restart_obj, domain, restart_frame)
+ restart_frame = restart_frame + 1
+ end if
+ endif
end do
end subroutine mpas_core_run
@@ -193,21 +241,210 @@
call timestep(domain, dt)
- if (mod(itimestep, config_stats_interval) == 0) then
- block_ptr => domain % blocklist
- if(associated(block_ptr % next)) then
- write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- 'that there is only one block per processor.'
- end if
+ if (config_stats_interval > 0) then
+ if (mod(itimestep, config_stats_interval) == 0) then
+ block_ptr => domain % blocklist
+ if (associated(block_ptr % next)) then
+ write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+ 'that there is only one block per processor.'
+ end if
- call timer_start("global diagnostics")
- call computeGlobalDiagnostics(domain % dminfo, &
- block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
- itimestep, dt)
- call timer_stop("global diagnostics")
+ call timer_start("global diagnostics")
+ call computeGlobalDiagnostics(domain % dminfo, &
+ block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global diagnostics")
+ end if
end if
end subroutine mpas_timestep
+
+
+subroutine init_ZLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell
+ real (kind=RKIND), dimension(:), pointer :: &
+ hZLevel, zMidZLevel, zTopZLevel, &
+ hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ hZLevel => block % mesh % hZLevel % array
+ zMidZLevel => block % mesh % zMidZLevel % array
+ zTopZLevel => block % mesh % zTopZLevel % array
+ nVertLevels = block % mesh % nVertLevels
+ hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
+ hRatioZLevelK => block % mesh % hRatioZLevelK % array
+ hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
+
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ ! Upon restart, the correct hZLevel should be in restart.nc
+ if (.not. config_do_restart) hZLevel = h(:,1)
+
+ ! hZLevel should be in the grid.nc and restart.nc file,
+ ! and h for k=1 must be specified there as well.
+
+ zTopZLevel(1) = 0.0
+ do k = 1,nVertLevels
+ zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
+ zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
+ end do
+
+ hMeanTopZLevel(1) = 0.0
+ hRatioZLevelK(1) = 0.0
+ hRatioZLevelKm1(1) = 0.0
+ do k = 2,nVertLevels
+ hMeanTopZLevel(k) = 0.5*(hZLevel(k-1) + hZLevel(k))
+ hRatioZLevelK(k) = 0.5*hZLevel(k)/hMeanTopZLevel(k)
+ hRatioZLevelKm1(k) = 0.5*hZLevel(k-1)/hMeanTopZLevel(k)
+ end do
+
+ block => block % next
+
+ end do
+
+end subroutine init_ZLevel
+
+
+subroutine compute_maxLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
+ real (kind=RKIND) :: centerx, centery
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+ integer, dimension(:), pointer :: &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell, &
+ boundaryVertex, verticesOnEdge
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+ maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+ maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ verticesOnEdge => block % mesh % verticesOnEdge % array
+ boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryCell => block % mesh % boundaryCell % array
+ boundaryVertex => block % mesh % boundaryVertex % array
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+ vertexDegree = block % mesh % vertexDegree
+
+ ! for z-grids, maxLevelCell should be in input state
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
+ maxLevelCell(1:nCells) = nVertLevels
+ endif
+ maxLevelCell(nCells+1) = 0
+
+ ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeTop(nEdges+1) = 0
+
+ ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeBot(iEdge) = &
+ max( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeBot(nEdges+1) = 0
+
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexBot(nVertices+1) = 0
+
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
+
+ ! set boundary edge
+ boundaryEdge=1
+ do iEdge=1,nEdges
+ boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
+ end do
+
+ !
+ ! Find cells and vertices that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
+ boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
+ endif
+ end do
+ end do
+
+ block => block % next
+ end do
+
+ ! Note: We do not update halos on maxLevel* variables. I want the
+ ! outside edge of a halo to be zero on each processor.
+
+end subroutine compute_maxLevel
subroutine mpas_core_finalize(domain)
Modified: branches/land_ice/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- branches/land_ice/mpas/src/core_ocean/module_test_cases.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_ocean/module_test_cases.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -25,15 +25,6 @@
type (block_type), pointer :: block_ptr
type (dm_info) :: dminfo
- integer :: iTracer
- real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
- integer :: nCells, nEdges, nVertices, nVertLevels
-
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -95,123 +86,6 @@
block_ptr => block_ptr % next
end do
- ! Initialize z-level grid variables from h, read in from input file.
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- h => block_ptr % state % time_levs(1) % state % h % array
- u => block_ptr % state % time_levs(1) % state % u % array
- rho => block_ptr % state % time_levs(1) % state % rho % array
- tracers => block_ptr % state % 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
- latCell => block_ptr % mesh % latCell % array
- lonCell => block_ptr % mesh % lonCell % array
-
- hZLevel => block_ptr % mesh % hZLevel % array
- zMidZLevel => block_ptr % mesh % zMidZLevel % array
- zTopZLevel => block_ptr % mesh % zTopZLevel % array
-
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
- pi=3.1415
- ! Tracer blob in Central Pacific, away from boundaries:
- !latCenter=pi/16; lonCenter=9./8.*pi
-
- ! Tracer blob in Central Pacific, near boundaries:
- latCenter=pi*2./16; lonCenter=13./16.*pi
-
- if (config_vert_grid_type.eq.'zlevel') then
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- hZLevel = h(:,1)
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- 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
- call dmpar_abort(dminfo)
- endif
-
- ! Set tracers, if not done in grid.nc file
- !tracers = 0.0
- centerx = 1.0e6
- centery = 1.0e6
- dist = 2.0e5
- do iCell = 1,nCells
- dist = sqrt( (latCell(iCell)-latCenter)**2 + (lonCell(iCell)-lonCenter)**2)
- do iLevel = 1,nVertLevels
- ! for 20 layer test
- ! tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 5.0 ! temperature
- ! tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- ! for x3, 25 layer test
- !tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 10.0 ! temperature
- !tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
- ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
-
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = block_ptr % mesh % lonCell % array(iCell)
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = block_ptr % mesh % latCell % array(iCell)
-
- ! place tracer blob with radius dist at (centerx,centery)
- !if ( sqrt( (xCell(iCell)-centerx)**2 &
- ! + (yCell(iCell)-centery)**2) &
- ! .lt.dist) then
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- !else
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 0.0
- !endif
-
- rho(iLevel,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) &
- + 7.6e-4*tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell))
-
- enddo
- enddo
-
- endif
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min u_src max u_src ', &
- ' min h max h ',&
- ' hZLevel zMidZlevel', &
- ' zTopZlevel'
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,1:nEdges)), maxval(u(iLevel,1:nEdges)), &
- minval(u_src(iLevel,1:nEdges)), maxval(u_src(iLevel,1:nEdges)), &
- minval(h(iLevel,1:nCells)), maxval(h(iLevel,1:nCells)), &
- hZLevel(iLevel),zMidZlevel(iLevel),zTopZlevel(iLevel)
- enddo
-
- print '(10a)', 'itracer ilevel min tracer max tracer'
- do iTracer=1,block_ptr % state % time_levs(1) % state % num_tracers
- do iLevel = 1,nVertLevels
- print '(2i5,20es12.4)', iTracer,ilevel, &
- minval(tracers(itracer,iLevel,1:nCells)), maxval(tracers(itracer,iLevel,1:nCells))
- enddo
- enddo
-
- block_ptr => block_ptr % next
- end do
-
-
end subroutine setup_sw_test_case
Modified: branches/land_ice/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_ocean/module_time_integration.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_ocean/module_time_integration.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -5,6 +5,7 @@
use constants
use dmpar
use vector_reconstruction
+ use spline_interpolation
contains
@@ -71,7 +72,7 @@
type (block_type), pointer :: block
type (state_type) :: provis
- integer :: rk_step
+ integer :: rk_step, iEdge
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
@@ -92,7 +93,7 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
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 % h % array(k,iCell)
end do
@@ -175,7 +176,7 @@
provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
provis % tracers % array(:,k,iCell) = ( &
block % state % time_levs(1) % state % h % array(k,iCell) * &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
@@ -204,7 +205,7 @@
+ rk_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
@@ -225,7 +226,7 @@
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
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 % h % array(k,iCell)
@@ -264,22 +265,23 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND) :: h_mom_eddy_visc2, h_mom_eddy_visc4
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- divergence, MontPot, pZLevel, zMidEdge, zTopEdge
+ MontPot, wTop, divergence
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
@@ -303,12 +305,10 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- ke_edge => s % ke_edge % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
- pZLevel => s % pZLevel % array
- zTopEdge => s % zTopEdge % array
- zMidEdge => s % zMidEdge % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -329,67 +329,75 @@
fEdge => grid % fEdge % array
zMidZLevel => grid % zMidZLevel % array
zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
tend_h => tend % h % array
tend_u => tend % u % array
nCells = grid % nCells
nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
u_src => grid % u_src % array
- h_mom_eddy_visc2 = config_h_mom_eddy_visc2
- h_mom_eddy_visc4 = config_h_mom_eddy_visc4
+ !
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_h = 0.0
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
- tend_h(:,:) = 0.0
- do iEdge=1,nEdges
+ !
+ ! for z-level, only compute height tendency for top layer.
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- 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 <= 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
- end do
- end if
- end do
- do iCell=1,nCells
- do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
- end do
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,min(1,maxLevelEdgeTop(iEdge))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
+ end do
+ do iCell=1,nCells
+ tend_h(1,iCell) = tend_h(1,iCell) / areaCell(iCell)
+ end do
+
+ endif ! config_vert_grid_type
+
!
! height tendency: vertical advection term -d/dz(hw)
!
+ ! Vertical advection computed for top layer of a z grid only.
if (config_vert_grid_type.eq.'zlevel') then
-
do iCell=1,nCells
-
tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
-
- ! This next loop is to verify that h for levels below the first
- ! remain constant. At a later time this could be replaced
- ! by just tend_h(2:nVertLevels,:) = 0.0, and then there is
- ! no need to compute the horizontal tend_h term for k=2:nVertLevels
- ! on a zlevel grid, above.
- do k=2,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) &
- - (wTop(k,iCell) - wTop(k+1,iCell))
- end do
-
end do
endif ! coordinate type
@@ -401,30 +409,30 @@
!
! velocity tendency: vertical advection term -w du/dz
!
- allocate(w_dudzTopEdge(nVertLevels+1))
- w_dudzTopEdge(1) = 0.0
- w_dudzTopEdge(nVertLevels+1) = 0.0
if (config_vert_grid_type.eq.'zlevel') then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ allocate(w_dudzTopEdge(nVertLevels+1))
+ w_dudzTopEdge(1) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- do k=2,nVertLevels
- ! Average w from cell center to edge
- wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+ do k=2,maxLevelEdgeTop(iEdge)
+ ! Average w from cell center to edge
+ wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
- ! compute dudz at vertical interface with first order derivative.
- w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
- / (zMidZLevel(k-1) - zMidZLevel(k))
- end do
+ ! compute dudz at vertical interface with first order derivative.
+ w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
+ / (zMidZLevel(k-1) - zMidZLevel(k))
+ end do
+ w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
- ! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,nVertLevels
- tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
- enddo
- enddo
+ ! Average w*du/dz from vertical interface to vertical middle of cell
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+ enddo
+ enddo
+ deallocate(w_dudzTopEdge)
endif
- deallocate(w_dudzTopEdge)
!
! velocity tendency: pressure gradient
@@ -434,7 +442,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
@@ -443,10 +451,10 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pZLevel(k,cell2) &
- - pZLevel(k,cell1) )/dcEdge(iEdge)
+ - rho0Inv*( pressure(k,cell2) &
+ - pressure(k,cell1) )/dcEdge(iEdge)
end do
enddo
endif
@@ -454,16 +462,16 @@
!
! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="black">abla^2 u
! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity )
- ! strictly only valid for h_mom_eddy_visc2 == constant
+ ! strictly only valid for config_h_mom_eddy_visc2 == constant
!
- if ( h_mom_eddy_visc2 > 0.0 ) then
+ if ( config_h_mom_eddy_visc2 > 0.0 ) then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
! is - </font>
<font color="gray">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
@@ -471,7 +479,7 @@
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = h_mom_eddy_visc2 * u_diffusion
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
@@ -483,9 +491,9 @@
! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="black">abla^4 u
! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
! applied recursively.
- ! strictly only valid for h_mom_eddy_visc4 == constant
+ ! strictly only valid for config_h_mom_eddy_visc4 == constant
!
- if ( h_mom_eddy_visc4 > 0.0 ) then
+ if ( config_h_mom_eddy_visc4 > 0.0 ) then
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_u(nVertLevels, nEdges+1))
@@ -501,12 +509,11 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
+ delsq_u(k,iEdge) = &
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
end do
end do
@@ -516,7 +523,7 @@
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- dcEdge(iEdge) * delsq_u(k,iEdge)
delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
@@ -525,7 +532,7 @@
end do
do iVertex=1,nVertices
r = 1.0 / areaTriangle(iVertex)
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
end do
end do
@@ -535,7 +542,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ delsq_u(k,iEdge)*dvEdge(iEdge)
delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
@@ -544,7 +551,7 @@
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
end do
end do
@@ -557,14 +564,14 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
u_diffusion = ( delsq_divergence(k,cell2) &
- delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
-( delsq_vorticity(k,vertex2) &
- delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - config_h_mom_eddy_visc4 * u_diffusion
end do
end do
@@ -582,7 +589,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
@@ -600,24 +607,34 @@
!
! velocity tendency: forcing and bottom drag
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! know the bottom edge with nonzero velocity and place the drag there.
+
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)
+ k = maxLevelEdgeTop(iEdge)
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
- tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- - 1.0e-3*u(nVertLevels,iEdge) &
- *sqrt(2.0*ke_edge(nVertLevels,iEdge))/h_edge(nVertLevels,iEdge)
+ ! efficiency note: it would be nice to avoid this
+ ! if within a do. This could be done with
+ ! k = max(maxLevelEdgeTop(iEdge),1)
+ ! and then tend_u(1,iEdge) is just not used for land cells.
- ! old bottom drag, just linear friction
- ! du/dt = u/tau where tau=100 days.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - u(nVertLevels,iEdge)/(100.0*86400.0)
+ if (k>0) then
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
+ endif
+
enddo
!
@@ -644,20 +661,23 @@
call dmpar_abort(dminfo)
endif
- allocate(fluxVertTop(1:nVertLevels+1))
+ allocate(fluxVertTop(nVertLevels+1))
fluxVertTop(1) = 0.0
- fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
- do k=2,nVertLevels
+
+ do k=2,maxLevelEdgeTop(iEdge)
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
- / (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
+ * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
enddo
- do k=1,nVertLevels
+ fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ (fluxVertTop(k) - fluxVertTop(k+1)) &
- /(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
+ / h_edge(k,iEdge)
enddo
+
end do
deallocate(fluxVertTop, vertViscTop)
@@ -680,39 +700,39 @@
type (mesh_type), intent(in) :: grid
integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
- nEdges, nCells, nVertLevels, num_tracers
+ nEdges, nCells, nCellsSolve, nVertLevels, num_tracers
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- real (kind=RKIND) :: flux, tracer_edge, r, dist
+ real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge, zMid, zTop
+ u,h,wTop, h_edge
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:,:), pointer :: boundaryEdge
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
+ hRatioZLevelK, hRatioZLevelKm1
+ real (kind=RKIND), dimension(:), allocatable:: vertDiffTop, tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
+ posZTopZLevel
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
- real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
-
u => s % u % array
h => s % h % array
boundaryCell=> grid % boundaryCell % array
wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
- zMid => s % zMid % array
- zTop => s % zTop % array
tend_tr => tend % tracers % array
@@ -721,10 +741,17 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
+ zMidZLevel => grid % zMidZLevel % array
+ hRatioZLevelK => grid % hRatioZLevelK % array
+ hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
boundaryEdge => grid % boundaryEdge % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
nEdges = grid % nEdges
nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
num_tracers = s % num_tracers
@@ -738,6 +765,11 @@
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
+ ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
+ ! tracer_edge at the boundary will also need to be defined for flux boundaries.
+
coef_3rd_order = 0.
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
@@ -747,16 +779,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,nVertLevels
- do iTracer=1,num_tracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
end do
- end if
+ end do
end do
else if (config_tracer_adv_order == 3) then
@@ -765,56 +795,52 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- 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) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,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) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(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) * u(k,iEdge) * h_edge(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. )
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(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) * u(k,iEdge) * h_edge(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. )
- end if
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
else if (config_tracer_adv_order == 4) then
@@ -823,46 +849,42 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- 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) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,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) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
endif ! if (config_tracer_adv_order == 2 )
@@ -871,44 +893,183 @@
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
- allocate(tracerTop(num_tracers,nVertLevels+1))
- tracerTop(:,1)=0.0
- tracerTop(:,nVertLevels+1)=0.0
- do iCell=1,grid % nCellsSolve
- if (config_vert_tracer_adv.eq.'centered') then
- do k=2,nVertLevels
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
- +tracers(iTracer,k ,iCell))/2.0
- end do
- end do
-
- elseif (config_vert_tracer_adv.eq.'upwind') then
- do k=2,nVertLevels
- if (wTop(k,iCell)>0.0) then
- upwindCell = k
- else
- upwindCell = k-1
- endif
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
+ if (config_vert_grid_type.eq.'zlevel') then
- endif
+ allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
- do k=1,nVertLevels
+ ! Tracers at the top and bottom boundary are assigned nearest
+ ! cell-centered value, regardless of tracer interpolation method.
+ ! wTop=0 at top and bottom sets the boundary condition.
+ do iCell=1,nCellsSolve
do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ tracerTop(iTracer,1,iCell) = tracers(iTracer,1,iCell)
+ tracerTop(iTracer,maxLevelCell(iCell)+1,iCell) = &
+ tracers(iTracer,maxLevelCell(iCell),iCell)
end do
end do
- enddo ! iCell
- deallocate(tracerTop)
+ if (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.2) then
+ ! Compute tracerTop using centered stencil, a simple average.
+
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ ( tracers(iTracer,k-1,iCell) &
+ +tracers(iTracer,k ,iCell))/2.0
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.3) then
+
+ ! Compute tracerTop using 3rd order stencil. This is the same
+ ! as 4th order, but includes upwinding.
+
+ ! Hardwire flux3Coeff at 1.0 for now. Could add this to the
+ ! namelist, if desired.
+ flux3Coef = 1.0
+ do iCell=1,nCellsSolve
+ k=2
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ do k=3,maxLevelCell(iCell)-1
+ cSignWTop = sign(flux3Coef,wTop(k,iCell))
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ ( (-1.+ cSignWTop)*tracers(iTracer,k-2,iCell) &
+ +( 7.-3.*cSignWTop)*tracers(iTracer,k-1,iCell) &
+ +( 7.+3.*cSignWTop)*tracers(iTracer,k ,iCell) &
+ +(-1.- cSignWTop)*tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
+ end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.4) then
+
+ ! Compute tracerTop using 4rd order stencil [-1 7 7 -1]
+
+ do iCell=1,nCellsSolve
+ k=2
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ do k=3,maxLevelCell(iCell)-1
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ (- tracers(iTracer,k-2,iCell) &
+ +7.*tracers(iTracer,k-1,iCell) &
+ +7.*tracers(iTracer,k ,iCell) &
+ - tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
+ end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.2) then
+
+ ! Compute tracerTop using linear interpolation.
+
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! Note hRatio on the k side is multiplied by tracer at k-1
+ ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.3) then
+
+ ! Compute tracerTop using cubic spline interpolation.
+
+ allocate(tracer2ndDer(nVertLevels))
+ allocate(tracersIn(nVertLevels),tracersOut(nVertLevels), &
+ posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels-1))
+
+ ! For the ocean, zlevel coordinates are negative and decreasing,
+ ! but spline functions assume increasing, so flip to positive.
+
+ posZMidZLevel = -zMidZLevel(1:nVertLevels)
+ posZTopZLevel = -zTopZLevel(2:nVertLevels)
+
+ do iCell=1,nCellsSolve
+ ! mrp 110201 efficiency note: push tracer loop down
+ ! into spline subroutines to improve efficiency
+ do iTracer=1,num_tracers
+
+ ! Place data in arrays to avoid creating new temporary arrays for every
+ ! subroutine call.
+ tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
+
+ call CubicSplineCoefficients(posZMidZLevel, &
+ tracersIn, maxLevelCell(iCell), tracer2ndDer)
+
+ call InterpolateCubicSpline( &
+ posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
+ posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
+
+ tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
+
+ end do
+ end do
+
+ deallocate(tracer2ndDer)
+ deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
+
+ else
+
+ print *, 'Abort: Incorrect combination of ', &
+ 'config_vert_tracer_adv and config_vert_tracer_adv_order.'
+ print *, 'Use:'
+ print *, 'config_vert_tracer_adv=''stencil'' and config_vert_tracer_adv_order=2,3,4 or'
+ print *, 'config_vert_tracer_adv=''spline'' and config_vert_tracer_adv_order=2,3'
+ call dmpar_abort(dminfo)
+
+ endif ! vertical tracer advection method
+
+ do iCell=1,nCellsSolve
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
+ end do
+ end do
+ end do
+
+ deallocate(tracerTop)
+
+ endif ! ZLevel
+
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
@@ -927,7 +1088,7 @@
invAreaCell1 = 1.0/areaCell(cell1)
invAreaCell2 = 1.0/areaCell(cell2)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
tracer_turb_flux = config_h_tracer_eddy_diff2 &
@@ -970,7 +1131,7 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ dvEdge(iEdge)*h_edge(k,iEdge) &
@@ -985,9 +1146,9 @@
end do
- do iCell = 1, nCells
+ do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
end do
@@ -1001,21 +1162,20 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,grid % nVertLevels
- do iTracer=1,num_tracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 &
- *( delsq_tracer(iTracer,k,cell2) &
- - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * tracer_turb_flux
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 &
+ *( delsq_tracer(iTracer,k,cell2) &
+ - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * tracer_turb_flux
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
- - flux * invAreaCell1 * boundaryMask(k,iEdge)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
- + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
+ - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
+ + flux * invAreaCell2 * boundaryMask(k,iEdge)
- end do
+ enddo
enddo
-
end do
deallocate(delsq_tracer)
@@ -1048,29 +1208,30 @@
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
- fluxVertTop(:,nVertLevels+1) = 0.0
- do iCell=1,grid % nCellsSolve
- do k=2,nVertLevels
+ do iCell=1,nCellsSolve
+
+ do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
fluxVertTop(iTracer,k) = vertDiffTop(k) &
* (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
- / (zMid(k-1,iCell) -zMid(k,iCell))
+ * 2 / (h(k-1,iCell) + h(k,iCell))
enddo
enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
- do k=1,nVertLevels
- dist = zTop(k,iCell) - zTop(k+1,iCell)
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
+ ! This is h d/dz( fluxVertTop) but h and dz cancel, so
+ ! reduces to delta( fluxVertTop)
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + h(k,iCell)*(fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1))/dist
+ + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
enddo
enddo
enddo ! iCell loop
deallocate(fluxVertTop, vertDiffTop)
-
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
@@ -1102,25 +1263,30 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel, ssh
+ hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, ke_edge, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
+ rho, temperature, salinity
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -1144,15 +1310,8 @@
gradPVt => s % gradPVt % array
rho => s % rho % array
tracers => s % tracers % array
- zMid => s % zMid % array
- zTop => s % zTop % array
- zMidEdge => s % zMidEdge % array
- zTopEdge => s % zTopEdge % array
- p => s % p % array
- pZLevel => s % pZLevel % array
- pTop => s % pTop % array
MontPot => s % MontPot % array
- ssh => s % ssh % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1173,35 +1332,27 @@
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
deriv_two => grid % deriv_two % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
nCells = grid % nCells
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
!
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
-
-
- !
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
coef_3rd_order = 0.
if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
@@ -1209,141 +1360,129 @@
if (config_thickness_adv_order == 2) then
- do iEdge=1,grid % nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
end do
else if (config_thickness_adv_order == 3) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,grid % nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- 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) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(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) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(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
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(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
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
else if (config_thickness_adv_order == 4) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,grid % nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- 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) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(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) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
endif ! if(config_thickness_adv_order == 2)
!
- ! 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
+ ! set the velocity and height at dummy address
+ ! used -1e34 so error clearly occurs if these values are used.
!
- u(:,nEdges+1) = 0.0
+ u(:,nEdges+1) = -1e34
+ h(:,nCells+1) = -1e34
+ tracers(s % index_temperature,:,nCells+1) = -1e34
+ tracers(s % index_salinity,:,nCells+1) = -1e34
!
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- 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) <= nVertices) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
+ circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
+ end do
end do
do iVertex=1,nVertices
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
end do
end do
-
!
! Compute the divergence at each cell center
!
@@ -1351,259 +1490,213 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end if
+ do k=1,maxLevelEdgeBot(iEdge)
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
end do
do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
+ r = 1.0 / areaCell(iCell)
+ do k = 1,maxLevelCell(iCell)
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
enddo
!
! Compute kinetic energy in each cell
!
ke(:,:) = 0.0
- do iCell=1,nCells
- 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) * u(k,iEdge)**2.0
- end do
- end do
- do k=1,nVertLevels
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ enddo
+ end do
+ do iCell = 1,nCells
+ do k = 1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
- end do
+ enddo
+ enddo
!
- !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe <= nEdges) then
- do k = 1,nVertLevels
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end if
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k = 1,maxLevelEdgeTop(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
end do
end do
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
+ ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
+ ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ ke_edge = 0.0 !mrp remove 0 for efficiency
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
- end do
- else
- do k=1,nVertLevels
- ke_edge(k,iEdge) = 0.0
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
end do
!
! 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) > nCells) cycle VTX_LOOP
- end do
- do k=1,nVertLevels
+ do iVertex = 1,nVertices
+ do k=1,maxLevelVertexBot(iVertex)
h_vertex = 0.0
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
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
+ end do
-
!
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
+ pv_cell(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ do k = 1,maxLevelCell(iCell)
+ pv_cell(k,iCell) = pv_cell(k,iCell) &
+ + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
+ / areaCell(iCell)
+ enddo
enddo
enddo
!
! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
+ ! ( this computes pv_edge at all edges bounding real cells )
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- if(iEdge <= nEdges) then
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ do i=1,vertexDegree
+ iEdge = edgesOnVertex(i,iVertex)
+ do k=1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
- endif
end do
end do
!
- ! Modify PV edge with upstream bias.
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradPVn for all edges bounding real cells )
!
+ gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
+ - pv_cell(k,cellsOnEdge(1,iEdge))) &
+ / dcEdge(iEdge)
enddo
enddo
-
!
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
- pv_cell(:,:) = 0.0
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- 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
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
+ - pv_vertex(k,verticesOnEdge(1,iEdge))) &
+ /dvEdge(iEdge)
+ enddo
enddo
!
- ! 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) <= 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
-
-
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ do k = 1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) &
+ - 0.5 * dt* ( u(k,iEdge) * gradPVn(k,iEdge) &
+ + v(k,iEdge) * gradPVt(k,iEdge) )
enddo
enddo
!
- ! Compute sea surface height
+ ! equation of state
!
- do iCell=1,nCells
- ssh(iCell) = h(1,iCell) - hZLevel(1)
- enddo
+ ! For an isopycnal model, density should remain constant.
+ if (config_vert_grid_type.eq.'zlevel') then
+ call equation_of_state(s,grid)
+ endif
!
- ! equation of state
+ ! Pressure
+ ! This section must be after computing rho
!
- ! For a isopycnal model, density should remain constant.
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ ! For Isopycnal model.
+ ! Compute pressure at top of each layer, and then
+ ! Montgomery Potential.
+ allocate(pTop(nVertLevels))
do iCell=1,nCells
- do k=1,nVertLevels
- ! Linear equation of state, for the time being
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(s % index_temperature,k,iCell) &
- + 7.6e-4*tracers(s % index_salinity,k,iCell))
- end do
- end do
- endif
+ ! assume atmospheric pressure at the surface is zero for now.
+ pTop(1) = 0.0
+ ! For isopycnal mode, p is the Montgomery Potential.
+ ! At top layer it is g*SSH, where SSH may be off by a
+ ! constant (ie, h_s can be relative to top or bottom)
+ MontPot(1,iCell) = gravity &
+ * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
- ! For Isopycnal model.
- ! Compute mid- and bottom-depth of each layer, from bottom up.
- ! Then compute mid- and bottom-pressure of each layer, and
- ! Montgomery Potential, from top down
- !
- do iCell=1,nCells
+ do k=2,nVertLevels
+ pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
- ! 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.
- zTop(nVertLevels+1,iCell) = h_s(iCell)
- do k=nVertLevels,1,-1
- zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
- zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
- end do
+ ! from delta M = p delta / rho
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ end do
- ! assume atmospheric pressure at the surface is zero for now.
- pTop(1,iCell) = 0.0
- do k=1,nVertLevels
- delta_p = rho(k,iCell)*gravity* h(k,iCell)
- p(k ,iCell) = pTop(k,iCell) + 0.5*delta_p
- pTop(k+1,iCell) = pTop(k,iCell) + delta_p
- end do
+ end do
+ deallocate(pTop)
- MontPot(1,iCell) = gravity * zTop(1,iCell)
- do k=2,nVertLevels
- ! from delta M = p delta / rho
- MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pTop(k,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
+ elseif (config_vert_grid_type.eq.'zlevel') then
- end do
+ ! For z-level model.
+ ! Compute pressure at middle of each level.
+ ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
+ ! pressure at middle of layer including SSH.
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if(cell1<=nCells .and. cell2<=nCells) then
- do k=1,nVertLevels
- zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
- zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
- enddo
- zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
- + zTop(nVertLevels+1,cell2))/2.0
- endif
- enddo
+ do iCell=1,nCells
+ ! compute pressure for z-level coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
+ pressure(1,iCell) = rho(1,iCell)*gravity &
+ * (h(1,iCell)-0.5*hZLevel(1))
- ! For z-level model.
- ! Compute pressure at middle of each level.
- ! pZLevel and p should only differ at k=1, where p is
- ! pressure at middle of layer including SSH, and pZLevel is
- ! pressure at a depth of hZLevel(1)/2.
- !
- do iCell=1,nCells
- ! compute pressure for z-level coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pZLevel(1,iCell) = rho(1,iCell)*gravity &
- * (h(1,iCell)-0.5*hZLevel(1))
- do k=2,nVertLevels
- pZLevel(k,iCell) = pZLevel(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
- + rho(k ,iCell)*hZLevel(k ))
- end do
+ do k=2,maxLevelCell(iCell)
+ pressure(k,iCell) = pressure(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ + rho(k ,iCell)*hZLevel(k ))
+ end do
- end do
+ end do
- ! compute vertical velocity
+ endif
+
+ !
+ ! vertical velocity through layer interface
+ !
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
wTop=0.0
@@ -1614,43 +1707,32 @@
! Compute div(u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- allocate(div_u(nVertLevels,nCells))
+ allocate(div_u(nVertLevels,nCells+1))
div_u(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell1) = div_u(k,cell1) + flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell2) = div_u(k,cell2) - flux
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,maxLevelEdgeBot(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge)
+ div_u(k,cell1) = div_u(k,cell1) + flux
+ div_u(k,cell2) = div_u(k,cell2) - flux
+ end do
end do
do iCell=1,nCells
- do k=1,nVertLevels
- div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) &
+ - div_u(k,iCell)/areaCell(iCell)*h(k,iCell)
end do
-
- ! Vertical velocity at bottom is zero.
- ! this next line can be set permanently somewhere upon startup.
- wTop(nVertLevels+1,iCell) = 0.0
- do k=nVertLevels,1,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
- end do
-
end do
deallocate(div_u)
endif
-
end subroutine compute_solve_diagnostics
@@ -1696,4 +1778,261 @@
end subroutine enforce_boundaryEdge
+
+ subroutine equation_of_state(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ if (config_eos_type.eq.'linear') then
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ elseif (config_eos_type.eq.'jm') then
+
+ call equation_of_state_jm(s, grid)
+
+ else
+ print *, ' Incorrect choice of config_eos_type:',&
+ config_eos_type
+ call dmpar_abort(dminfo)
+ endif
+
+ end subroutine equation_of_state
+
+
+ subroutine equation_of_state_jm(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND) :: p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ p = pRefEOS(k)
+ p2 = pRefEOS(k)*pRefEOS(k)
+
+ SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p)
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p)
+
+ rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS)
+
+ end subroutine equation_of_state_jm
+
end module time_integration
Modified: branches/land_ice/mpas/src/core_sw/Registry
===================================================================
--- branches/land_ice/mpas/src/core_sw/Registry 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_sw/Registry 2011-05-10 20:41:41 UTC (rev 826)
@@ -15,10 +15,13 @@
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
+namelist logical sw_model config_wind_stress false
+namelist logical sw_model config_bottom_drag false
namelist real sw_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
+namelist character io config_input_name grid.nc
+namelist character io config_output_name output.nc
+namelist character io config_restart_name restart.nc
+namelist character io config_decomp_file_prefix graph.info.part.
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -111,6 +114,7 @@
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 - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
Modified: branches/land_ice/mpas/src/core_sw/module_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_sw/module_time_integration.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/core_sw/module_time_integration.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -269,8 +269,10 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+ real (kind=RKIND) :: ke_edge
-
h => s % h % array
u => s % u % array
v => s % v % array
@@ -308,26 +310,20 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ u_src => grid % u_src % array
!
! Compute height tendency for each cell
!
tend_h(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- 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 <= 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
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
@@ -504,6 +500,28 @@
end if
+ ! Compute u (velocity) tendency from wind stress (u_src)
+ 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)
+ end do
+ endif
+
+ if (config_bottom_drag) then
+ do iEdge=1,grid % nEdges
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+ ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
+ + ke(1,cellsOnEdge(2,iEdge)))
+
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ - 1.0e-3*u(1,iEdge) &
+ *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+ end do
+ endif
+
end subroutine compute_tend
@@ -1092,10 +1110,7 @@
! 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) > nCells) cycle VTX_LOOP
- end do
+ do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex(k,iVertex) = 0.0
do i=1,grid % vertexDegree
@@ -1105,7 +1120,7 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
end do
- end do VTX_LOOP
+ end do
!
Modified: branches/land_ice/mpas/src/framework/Makefile
===================================================================
--- branches/land_ice/mpas/src/framework/Makefile 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/framework/Makefile 2011-05-10 20:41:41 UTC (rev 826)
@@ -31,7 +31,7 @@
module_dmpar.o: module_sort.o streams.o
-module_block_decomp.o: module_grid_types.o module_hash.o
+module_block_decomp.o: module_grid_types.o module_hash.o module_configure.o
module_io_input.o: module_grid_types.o module_dmpar.o module_block_decomp.o module_sort.o module_configure.o $(ZOLTANOBJ)
Modified: branches/land_ice/mpas/src/framework/module_block_decomp.F
===================================================================
--- branches/land_ice/mpas/src/framework/module_block_decomp.F 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/framework/module_block_decomp.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -18,6 +18,8 @@
subroutine block_decomp_cells_for_proc(dminfo, partial_global_graph_info, local_cell_list)
+ use configure
+
implicit none
type (dm_info), intent(in) :: dminfo
@@ -41,15 +43,15 @@
iunit = 50 + dminfo % my_proc_id
if (dminfo % nprocs < 10) then
- write(filename,'(a,i1)') 'graph.info.part.', dminfo % nprocs
+ write(filename,'(a,i1)') trim(config_decomp_file_prefix), dminfo % nprocs
else if (dminfo % nprocs < 100) then
- write(filename,'(a,i2)') 'graph.info.part.', dminfo % nprocs
+ write(filename,'(a,i2)') trim(config_decomp_file_prefix), dminfo % nprocs
else if (dminfo % nprocs < 1000) then
- write(filename,'(a,i3)') 'graph.info.part.', dminfo % nprocs
+ write(filename,'(a,i3)') trim(config_decomp_file_prefix), dminfo % nprocs
else if (dminfo % nprocs < 10000) then
- write(filename,'(a,i4)') 'graph.info.part.', dminfo % nprocs
+ write(filename,'(a,i4)') trim(config_decomp_file_prefix), dminfo % nprocs
else if (dminfo % nprocs < 100000) then
- write(filename,'(a,i5)') 'graph.info.part.', dminfo % nprocs
+ write(filename,'(a,i5)') trim(config_decomp_file_prefix), dminfo % nprocs
end if
open(unit=iunit, file=trim(filename), form='formatted', status='old', iostat=istatus)
Modified: branches/land_ice/mpas/src/operators/Makefile
===================================================================
--- branches/land_ice/mpas/src/operators/Makefile 2011-05-10 18:38:01 UTC (rev 825)
+++ branches/land_ice/mpas/src/operators/Makefile 2011-05-10 20:41:41 UTC (rev 826)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_vector_reconstruction.o
+OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o
all: operators
@@ -9,6 +9,7 @@
module_vector_reconstruction.o: module_RBF_interpolation.o
module_RBF_interpolation.o:
+module_spline_interpolation:
clean:
$(RM) *.o *.mod *.f90 libops.a
Copied: branches/land_ice/mpas/src/operators/module_spline_interpolation.F (from rev 825, trunk/mpas/src/operators/module_spline_interpolation.F)
===================================================================
--- branches/land_ice/mpas/src/operators/module_spline_interpolation.F (rev 0)
+++ branches/land_ice/mpas/src/operators/module_spline_interpolation.F 2011-05-10 20:41:41 UTC (rev 826)
@@ -0,0 +1,427 @@
+module spline_interpolation
+
+ implicit none
+
+ private
+
+ public :: CubicSplineCoefficients, InterpolateCubicSpline, &
+ IntegrateCubicSpline, IntegrateColumnCubicSpline, InterpolateLinear, &
+ TestInterpolate
+
+! Short Descriptions:
+
+! CubicSplineCoefficients: Compute second derivatives at nodes.
+! This must be run before any of the other cubic spine functions.
+
+! InterpolateCubicSpline: Compute cubic spline interpolation.
+
+! IntegrateCubicSpline: Compute a single integral from spline data.
+
+! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
+
+! InterpolateLinear: Compute linear interpolation.
+
+! TestInterpolate: Test spline interpolation subroutines.
+
+ contains
+
+ subroutine CubicSplineCoefficients(x,y,n,y2ndDer)
+
+! Given arrays x(1:n) and y(1:n) containing a function,
+! i.e., y(i) = f(x(i)), with x monotonically increasing
+! this routine returns an array y2ndDer(1:n) that contains
+! the second derivatives of the interpolating function at x(1:n).
+! This routine uses boundary conditions for a natural spline,
+! with zero second derivative on that boundary.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y ! value at nodes
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out), dimension(n) :: &
+ y2ndDer ! dy^2/dx^2 at each node
+
+! local variables:
+
+ integer :: i
+ real(kind=RKIND) :: &
+ temp,xRatio,a(n)
+
+ y2ndDer(1)=0.0
+ y2ndDer(n)=0.0
+ a(1)=0.0
+
+ do i=2,n-1
+ xRatio=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ temp=1.0/(2.0+xRatio*y2ndDer(i-1))
+ y2ndDer(i)=temp*(xRatio-1.0)
+ a(i) = temp*(6.0*((y(i+1)-y(i))/(x(i+1)-x(i)) &
+ -(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1)) &
+ -xRatio*a(i-1))
+ enddo
+
+ do i=n-1,1,-1
+ y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
+ enddo
+
+ end subroutine CubicSplineCoefficients
+
+
+ subroutine InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns the
+! cubic-spline interpolated values of yOut(1:nOut) at xOut(1:nOut).
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+! last values of x.
+
+! INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y, &! interpolation variable, input grid
+ y2ndDer ! 2nd derivative of y at nodes
+
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
+
+ integer, intent(in) :: &
+ n, &! number of nodes, input grid
+ nOut ! number of nodes, output grid
+
+! OUTPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
+
+! local variables:
+
+ integer :: &
+ kIn, kOut ! counters
+
+ real (kind=RKIND) :: &
+ a, b, h
+
+ kOut = 1
+
+ kInLoop: do kIn = 1,n-1
+
+ h = x(kIn+1)-x(kIn)
+
+ do while(xOut(kOut) < x(kIn+1))
+
+ a = (x(kIn+1)-xOut(kOut))/h
+ b = (xOut(kOut)-x (kIn) )/h
+ yOut(kOut) = a*y(kIn) + b*y(kIn+1) &
+ + ((a**3-a)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
+ *(h**2)/6.0
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+end subroutine InterpolateCubicSpline
+
+
+subroutine IntegrateCubicSpline(x,y,y2ndDer,n,x1,x2,y_integral)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns y_integral,
+! the integral of y from x1 to x2. The integration formula was
+! created by analytically integrating a cubic spline between each node.
+! This subroutine assumes that x is monotonically increasing, and
+! that x1 < x2.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y, &! value at nodes
+ y2ndDer ! dy^2/dx^2 at each node
+ real(kind=RKIND), intent(in) :: &
+ x1,x2 ! limits of integration
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out) :: &
+ y_integral ! integral of y
+
+! local variables:
+
+ integer :: i,j,k
+ real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
+
+ if (x1<x(1).or.x2>x(n).or.x1>x2) then
+ print *, 'error on integration bounds'
+ endif
+
+ y_integral = 0.0
+ eps1 = 1e-14*x2
+
+ do j=1,n-1 ! loop through sections
+ ! section x(j) ... x(j+1)
+
+ if (x2<=x(j) +eps1) exit
+ if (x1>=x(j+1)-eps1) cycle
+
+ h = x(j+1) - x(j)
+ h2 = h**2
+
+ ! left side:
+ if (x1<x(j)) then
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ else
+ A2 = (x(j+1)-x1 )**2/h2
+ B2 = (x1 -x(j))**2/h2
+ F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ ! right side:
+ if (x2>x(j+1)) then
+ F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
+ else
+ A2 = (x(j+1)-x2 )**2/h2
+ B2 = (x2 -x(j))**2/h2
+ F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ y_integral = y_integral + F2 - F1
+
+ enddo ! j
+
+ end subroutine IntegrateCubicSpline
+
+
+ subroutine IntegrateColumnCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,y_integral, nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns
+! y_integral(1:nOut), the integral of y.
+! This is a cumulative integration, so that
+! y_integral(j) holds the integral of y from x(1) to xOut(j).
+! The integration formula was created by analytically integrating a
+! cubic spline between each node.
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n, &! number of nodes
+ nOut ! number of output locations to compute integral
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y, &! value at nodes
+ y2ndDer ! dy^2/dx^2 at each node
+ real(kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! output locations to compute integral
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), dimension(nOut), intent(out) :: &
+ y_integral ! integral from 0 to xOut
+
+! local variables:
+
+ integer :: i,j,k
+ real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
+
+ y_integral = 0.0
+ j = 1
+ h = x(j+1) - x(j)
+ h2 = h**2
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
+
+ k_loop: do k = 1,nOut
+
+ if (k>1) y_integral(k) = y_integral(k-1)
+
+ do while(xOut(k) > x(j+1)-eps1)
+ F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
+
+ y_integral(k) = y_integral(k) + F2 - F1
+ j = j+1
+ h = x(j+1) - x(j)
+ h2 = h**2
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ if (abs(xOut(k) - x(j+1))<eps1) cycle k_loop
+ enddo
+
+ A2 = (x(j+1) - xOut(k))**2/h2
+ B2 = (xOut(k) - x(j) )**2/h2
+ F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+
+ y_integral(k) = y_integral(k) + F2 - F1
+
+ if (k < nOut) then
+ A2 = (x(j+1) -xOut(k))**2/h2
+ B2 = (xOut(k) -x(j) )**2/h2
+ F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ enddo k_loop
+
+ end subroutine IntegrateColumnCubicSpline
+
+
+ subroutine InterpolateLinear( &
+ x,y,n, &
+ xOut,yOut,nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! this routine returns the linear interpolated values of yOut(1:nOut)
+! at xOut(1:nOut).
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+! last values of x.
+
+! !INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y ! interpolation variable, input grid
+
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
+
+ integer, intent(in) :: &
+ N, &! number of nodes, input grid
+ NOut ! number of nodes, output grid
+
+! !OUTPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
+
+!-----------------------------------------------------------------------
+!
+! local variables
+!
+!-----------------------------------------------------------------------
+
+ integer :: &
+ kIn, kOut ! counters
+
+ kOut = 1
+
+ kInLoop: do kIn = 1,n-1
+
+ do while(xOut(kOut) < x(kIn+1))
+
+ yOut(kOut) = y(kIn) &
+ + (y(kIn+1)-y(kIn)) &
+ /(x(kIn+1) -x(kIn) ) &
+ *(xOut(kOut) -x(kIn) )
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+ end subroutine InterpolateLinear
+
+
+ subroutine TestInterpolate
+
+! Test function to show how to operate the cubic spline subroutines
+
+ integer, parameter :: &
+ n = 10
+ real (kind=RKIND), dimension(n) :: &
+ y, x, y2ndDer
+
+ integer, parameter :: &
+ nOut = 100
+ real (kind=RKIND), dimension(nOut) :: &
+ yOut, xOut
+
+ integer :: &
+ k
+
+!-----------------------------------------------------------------------
+!
+! Create x, y, xOut
+!
+!-----------------------------------------------------------------------
+
+ do k=1,n
+ x(k) = k-4
+ ! trig function:
+ y(k) = sin(x(k)/2)
+ enddo
+
+ do k=1,nOut
+ xOut(k) = x(1) + k/(nOut+1.0)*(x(n)-x(1))
+ enddo
+
+!-----------------------------------------------------------------------
+!
+! Interpolate
+!
+!-----------------------------------------------------------------------
+
+ ! First, compute second derivative values at each node, y2ndDer.
+ call CubicSplineCoefficients(x,y,n,y2ndDer)
+
+ ! Compute interpolated values yOut.
+ call InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,1)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = [',y,'];'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+
+ ! Compute interpolated values yOut.
+ call IntegrateColumnCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,2)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = 2*cos(-3/2) -2*cos(x/2);'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+
+ end subroutine TestInterpolate
+
+end module spline_interpolation
+
</font>
</pre>