<p><b>ringler@lanl.gov</b> 2011-05-28 18:18:43 -0600 (Sat, 28 May 2011)</p><p><br>
add the same mesh scaling code as in the truck shallow-water model.<br>
<br>
this allows the del2/del4 mixing of u and scalars to vary with grid spacing.<br>
<br>
the namelist parameter (config_h_ScaleWithMesh) has a default value of false.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-28 23:34:25 UTC (rev 858)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-29 00:18:43 UTC (rev 859)
@@ -22,6 +22,7 @@
namelist integer timestep_higdon config_n_bcl_iter_end 4
namelist integer timestep_higdon config_n_btr_subcycles 10
namelist logical timestep_higdon config_compute_tr_midstage true
+namelist logical sw_model config_h_ScaleWithMesh false
namelist real hmix config_h_mom_eddy_visc2 0.0
namelist real hmix config_h_mom_eddy_visc4 0.0
namelist real hmix config_h_tracer_eddy_diff2 0.0
@@ -96,6 +97,10 @@
var persistent real zVertex ( nVertices ) 0 iro zVertex mesh - -
var persistent integer indexToVertexID ( nVertices ) 0 iro indexToVertexID mesh - -
+var persistent real meshDensity ( nCells ) 0 iro meshDensity mesh - -
+var persistent real meshScalingDel2 ( nEdges ) 0 ro meshScalingDel2 mesh - -
+var persistent real meshScalingDel4 ( nEdges ) 0 ro meshScalingDel4 mesh - -
+
var persistent integer cellsOnEdge ( TWO nEdges ) 0 iro cellsOnEdge mesh - -
var persistent integer nEdgesOnCell ( nCells ) 0 iro nEdgesOnCell mesh - -
var persistent integer nEdgesOnEdge ( nEdges ) 0 iro nEdgesOnEdge mesh - -
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-05-28 23:34:25 UTC (rev 858)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-05-29 00:18:43 UTC (rev 859)
@@ -80,6 +80,7 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+ call compute_mesh_scaling(mesh)
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
@@ -467,4 +468,37 @@
end subroutine mpas_core_finalize
+
+ subroutine compute_mesh_scaling(mesh)
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer :: iEdge, cell1, cell2
+ real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
+
+ meshDensity => mesh % meshDensity % array
+ meshScalingDel2 => mesh % meshScalingDel2 % array
+ meshScalingDel4 => mesh % meshScalingDel4 % array
+
+ !
+ ! Compute the scaling factors to be used in the del2 and del4 dissipation
+ !
+ meshScalingDel2(:) = 1.0
+ meshScalingDel4(:) = 1.0
+ if (config_h_ScaleWithMesh) then
+ do iEdge=1,mesh%nEdges
+ cell1 = mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = mesh % cellsOnEdge % array(2,iEdge)
+ meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
+ meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
+ end do
+ end if
+
+ end subroutine compute_mesh_scaling
+
end module mpas_core
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-28 23:34:25 UTC (rev 858)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-29 00:18:43 UTC (rev 859)
@@ -1234,7 +1234,7 @@
upstream_bias, wTopEdge, rho0Inv, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zMidZLevel, zTopZLevel
+ zMidZLevel, zTopZLevel, meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
@@ -1306,6 +1306,9 @@
u_src => grid % u_src % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
!
! velocity tendency: start accumulating tendency terms
!
@@ -1408,7 +1411,7 @@
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
+ u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
@@ -1499,8 +1502,10 @@
- delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
-( delsq_vorticity(k,vertex2) &
- delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) - config_h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
end do
end do
@@ -1719,7 +1724,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
- hRatioZLevelK, hRatioZLevelKm1
+ hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
posZTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
@@ -1762,6 +1767,10 @@
nVertLevels = grid % nVertLevels
num_tracers = s % num_tracers
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+
deriv_two => grid % deriv_two % array
if(config_restoreTS) then
@@ -2105,7 +2114,7 @@
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 &
+ tracer_turb_flux = meshScalingDel2(iEdge) * config_h_tracer_eddy_diff2 &
*( tracers(iTracer,k,cell2) &
- tracers(iTracer,k,cell1))/dcEdge(iEdge)
@@ -2178,7 +2187,7 @@
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 &
+ tracer_turb_flux = meshScalingDel4(iEdge) * config_h_tracer_eddy_diff4 &
*( delsq_tracer(iTracer,k,cell2) &
- delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * tracer_turb_flux
</font>
</pre>