<p><b>qchen3@fsu.edu</b> 2010-11-11 18:40:13 -0700 (Thu, 11 Nov 2010)</p><p>TRUNK COMMIT<br>
<br>
I am merging the branch swmodel_del4 back to the trunk. The branch was created to implement the del2 and del4 diffusions for both the momentum and the tracers. The code has been evaluated by Todd Ringler and Mark Petersen. A summary of changes is attached.<br>
<br>
CHANGES<br>
1. In core_sw/module_time_integration.F, in the subroutine compute_tend, routines were added to implement the del2 and del4 diffusions. The original code inside the LANL_FORMULATION for del2 diffusion has been removed. <br>
<br>
2. In core_sw/module_time_integration.F, in the subroutine compute_scalar_tend, routines were added to implement the del2 and del4 diffusions for the tracers. <br>
<br>
3. In support for the previous two changes, the following new variables were added to core_sw/Registry and namelist.input.sw:<br>
config_h_mom_eddy_visc2<br>
config_h_mom_eddy_visc4<br>
config_h_tracer_eddy_diff2<br>
config_h_tracer_eddy_diff4<br>
<br>
4. Also in support for Changes 1 & 2, we enforced conditional halo update on divergence and vorticity, at the beginning of the rk4 subroutine in core_sw/module_time_integration.F. We are not absolutely sure whether these updates are necessary. My testings show that they do make a difference.<br>
<br>
5. I added the gfortran-serial entry in the main Makefile<br>
<br>
6. On line 277 of core_sw/module_global_dynamics.F, '1i' was changed to '1i0' so that gfortran can compile. It seems that my version of gfortran (v4.4) insists on a non-negative digit as the field width after the descriptor 'i'. <br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/Makefile
===================================================================
--- trunk/mpas/Makefile 2010-11-12 00:14:40 UTC (rev 608)
+++ trunk/mpas/Makefile 2010-11-12 01:40:13 UTC (rev 609)
@@ -104,6 +104,18 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+gfortran-serial:
+ ( make all \
+ "FC = gfortran" \
+ "CC = gcc" \
+ "SFC = gfortran" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
+ "CFLAGS = -O3 -m64" \
+ "LDFLAGS = -O3 -m64" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
g95:
( make all \
"FC = mpif90" \
Modified: trunk/mpas/namelist.input.sw
===================================================================
--- trunk/mpas/namelist.input.sw 2010-11-12 00:14:40 UTC (rev 608)
+++ trunk/mpas/namelist.input.sw 2010-11-12 01:40:13 UTC (rev 609)
@@ -5,7 +5,10 @@
config_ntimesteps = 7500
config_output_interval = 500
config_stats_interval = 0
- config_visc = 0.0
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
config_tracer_adv_order = 2
config_positive_definite = .false.
Modified: trunk/mpas/src/core_sw/Registry
===================================================================
--- trunk/mpas/src/core_sw/Registry 2010-11-12 00:14:40 UTC (rev 608)
+++ trunk/mpas/src/core_sw/Registry 2010-11-12 01:40:13 UTC (rev 609)
@@ -7,7 +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 real sw_model config_visc 0.0
+namelist real sw_model config_h_mom_eddy_visc2 0.0
+namelist real sw_model config_h_mom_eddy_visc4 0.0
+namelist real sw_model config_h_tracer_eddy_diff2 0.0
+namelist real sw_model config_h_tracer_eddy_diff4 0.0
namelist integer sw_model config_thickness_adv_order 2
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
@@ -141,3 +144,4 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
+
Modified: trunk/mpas/src/core_sw/module_global_diagnostics.F
===================================================================
--- trunk/mpas/src/core_sw/module_global_diagnostics.F 2010-11-12 00:14:40 UTC (rev 608)
+++ trunk/mpas/src/core_sw/module_global_diagnostics.F 2010-11-12 01:40:13 UTC (rev 609)
@@ -274,7 +274,7 @@
else
open(fileID, file='GlobalIntegrals.txt',POSITION='append')
endif
- write(fileID,'(1i, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
+ write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
globalKineticEnergy, globalPotentialEnergy
close(fileID)
Modified: trunk/mpas/src/core_sw/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_sw/module_time_integration.F 2010-11-12 00:14:40 UTC (rev 608)
+++ trunk/mpas/src/core_sw/module_time_integration.F 2010-11-12 01:40:13 UTC (rev 609)
@@ -120,6 +120,16 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ end if
+
block => block % next
end do
@@ -253,11 +263,14 @@
circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: u_diffusion, visc
+ real (kind=RKIND) :: r, u_diffusion
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
- visc = config_visc
+
h => s % h % array
u => s % u % array
v => s % v % array
@@ -334,15 +347,6 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
-
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
- ! only valid for visc == constant
- !
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = visc * u_diffusion
-
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
@@ -352,12 +356,13 @@
tend_u(k,iEdge) = &
q &
- + u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) + &
gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
) / dcEdge(iEdge)
end do
end do
+
+
#endif
#ifdef NCAR_FORMULATION
@@ -386,6 +391,119 @@
end do
#endif
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for visc == constant
+ 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
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+ end if
+
+ !
+ ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="blue">abla^4 u
+ ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ ! applied recursively.
+ ! strictly only valid for h_mom_eddy_visc4 == constant
+ !
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
+
+ delsq_u(:,:) = 0.0
+
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ end do
+ end do
+
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end do
+ do iVertex=1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k=1,nVertLevels
+ delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ end do
+ end do
+
+ ! Divergence using </font>
<font color="blue">abla^2 u
+ delsq_divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ + delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+ - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ end do
+ end do
+
+ ! Compute - \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
+ 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
+
+ u_diffusion = ( delsq_divergence(k,cell2) &
+ - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( delsq_vorticity(k,vertex2) &
+ - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = config_h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
+ end if
+
end subroutine compute_tend
@@ -405,7 +523,12 @@
type (mesh_type), intent(in) :: grid
integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
+
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
@@ -422,6 +545,7 @@
tracers => s % tracers % array
cellsOnEdge => grid % cellsOnEdge % array
boundaryCell=> grid % boundaryCell % array
+ boundaryEdge=> grid % boundaryEdge % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
@@ -429,6 +553,7 @@
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
+
tracer_tend(:,:,:) = 0.0
if (config_tracer_adv_order == 2) then
@@ -556,6 +681,108 @@
endif ! if (config_tracer_adv_order == 2 )
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="blue">abla \phi)
+ !
+ if ( config_h_tracer_eddy_diff2 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1.0/areaCell(cell1)
+ invAreaCell2 = 1.0/areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+ tracer_turb_flux = config_h_tracer_eddy_diff2 &
+ *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
+
+ ! div(h \kappa_2 </font>
<font color="blue">abla \phi) at cell center
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
+ end do
+ end do
+
+ end do
+
+ deallocate(boundaryMask)
+
+ end if
+
+ !
+ ! tracer tendency: del4 horizontal tracer diffusion, &
+ ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="blue">abla \phi)])
+ !
+ if ( config_h_tracer_eddy_diff4 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ ! first del2: div(h </font>
<font color="blue">abla \phi) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
+ - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ end do
+ end do
+
+ end do
+
+ do iCell = 1, grid % nCells
+ r = 1.0 / grid % areaCell % array(iCell)
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
+ end do
+ end do
+
+ ! second del2: div(h </font>
<font color="blue">abla [delsq_tracer]) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
+ invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ 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
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ end do
+ enddo
+
+ end do
+
+ deallocate(delsq_tracer)
+ deallocate(boundaryMask)
+
+ end if
+
end subroutine compute_scalar_tend
</font>
</pre>