<p><b>qchen3@fsu.edu</b> 2010-09-17 11:59:55 -0600 (Fri, 17 Sep 2010)</p><p></p><hr noshade><pre><font color="gray">Modified: branches/swmodel_del4/Makefile
===================================================================
--- branches/swmodel_del4/Makefile 2010-09-17 16:40:45 UTC (rev 505)
+++ branches/swmodel_del4/Makefile 2010-09-17 17:59:55 UTC (rev 506)
@@ -96,18 +96,6 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
-ifort-serial:
- ( make all \
- "FC = ifort" \
- "CC = gcc" \
- "SFC = ifort" \
- "SCC = gcc" \
- "FFLAGS = -real-size 64 -O3" \
- "CFLAGS = -O3 -m64" \
- "LDFLAGS = -O3" \
- "CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
-
gfortran:
( make all \
"FC = mpif90" \
@@ -159,7 +147,7 @@
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf -lcurl
+LIBS = -L$(NETCDF)/lib -lnetcdf
RM = rm -f
CPP = cpp -C -P -traditional
Modified: branches/swmodel_del4/namelist.input.sw
===================================================================
--- branches/swmodel_del4/namelist.input.sw 2010-09-17 16:40:45 UTC (rev 505)
+++ branches/swmodel_del4/namelist.input.sw 2010-09-17 17:59:55 UTC (rev 506)
@@ -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 = 60000000000000.0
+ config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff4 = 0.0
/
&io
Modified: branches/swmodel_del4/src/core_sw/Registry
===================================================================
--- branches/swmodel_del4/src/core_sw/Registry 2010-09-17 16:40:45 UTC (rev 505)
+++ branches/swmodel_del4/src/core_sw/Registry 2010-09-17 17:59:55 UTC (rev 506)
@@ -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 character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
Modified: branches/swmodel_del4/src/core_sw/module_global_diagnostics.F
===================================================================
--- branches/swmodel_del4/src/core_sw/module_global_diagnostics.F 2010-09-17 16:40:45 UTC (rev 505)
+++ branches/swmodel_del4/src/core_sw/module_global_diagnostics.F 2010-09-17 17:59:55 UTC (rev 506)
@@ -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: branches/swmodel_del4/src/core_sw/module_time_integration.F
===================================================================
--- branches/swmodel_del4/src/core_sw/module_time_integration.F 2010-09-17 16:40:45 UTC (rev 505)
+++ branches/swmodel_del4/src/core_sw/module_time_integration.F 2010-09-17 17:59:55 UTC (rev 506)
@@ -247,11 +247,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, h_mom_eddy_visc2, h_mom_eddy_visc4
+ 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
@@ -289,7 +292,10 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ h_mom_eddy_visc2 = config_h_mom_eddy_visc2
+ h_mom_eddy_visc4 = config_h_mom_eddy_visc4
+
!
! Compute height tendency for each cell
!
@@ -328,15 +334,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)
@@ -346,12 +343,125 @@
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
+
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for visc == constant
+ if (h_mom_eddy_visc2 > 0.0) then
+ 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
+ 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
+ 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 (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 = 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
+
#endif
#ifdef NCAR_FORMULATION
@@ -399,8 +509,17 @@
type (grid_meta), intent(in) :: grid
integer :: iCell, iEdge, k, iTracer, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND) :: h_tracer_eddy_diff2, h_tracer_eddy_diff4, invAreaCell1, invAreaCell2, tracer_turb_flux
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
+
+ boundaryEdge => grid % boundaryEdge % array
+ h_tracer_eddy_diff2 = config_h_tracer_eddy_diff2
+ h_tracer_eddy_diff4 = config_h_tracer_eddy_diff4
+
tend % tracers % array(:,:,:) = 0.0
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
@@ -425,6 +544,119 @@
end do
end do
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="blue">abla \phi)
+ !
+ if ( 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 = 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
+ ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+ tracer_turb_flux = h_tracer_eddy_diff2 &
+ *( s % tracers % array(iTracer,k,cell2) &
+ - s % tracers % array(iTracer,k,cell1))/grid % dcEdge % array(iEdge)
+
+ ! div(h \kappa_2 </font>
<font color="blue">abla \phi) at cell center
+ flux = grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) &
+ * tracer_turb_flux * boundaryMask(k, iEdge)
+ tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) + flux * invAreaCell1
+ tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(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 ( 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 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ + grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) &
+ *(s % tracers % array(iTracer,k,cell2) - s % tracers % array(iTracer,k,cell1)) &
+ /grid % dcEdge % array(iEdge) * boundaryMask(k,iEdge)
+ delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
+ - grid % dvEdge % array(iEdge)*s % h_edge % array(k,iEdge) &
+ *(s % tracers % array(iTracer,k,cell2) - s % tracers % array(iTracer,k,cell1)) &
+ /grid % dcEdge % array(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="gray">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 = h_tracer_eddy_diff4 &
+ *( delsq_tracer(iTracer,k,cell2) &
+ - delsq_tracer(iTracer,k,cell1))/grid % dcEdge % array(iEdge)
+ flux = grid % dvEdge %array(iEdge) * tracer_turb_flux
+
+ tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) &
+ - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) &
+ + flux * invAreaCell2 * boundaryMask(k,iEdge)
+
+ end do
+ enddo
+
+ end do
+
+ deallocate(delsq_tracer)
+
+ end if
+
end subroutine compute_scalar_tend
@@ -665,11 +897,11 @@
!
! Modify PV edge with upstream bias.
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
- enddo
- enddo
+ !do iEdge = 1,nEdges
+ ! do k = 1,nVertLevels
+ ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ ! enddo
+ !enddo
!
@@ -707,11 +939,11 @@
! Modify PV edge with upstream bias.
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
- enddo
- enddo
+ !do iEdge = 1,nEdges
+ ! do k = 1,nVertLevels
+ ! pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ ! enddo
+ !enddo
!
! set pv_edge = fEdge / h_edge at boundary points
</font>
</pre>