<p><b>ringler@lanl.gov</b> 2011-05-25 15:04:15 -0600 (Wed, 25 May 2011)</p><p><br>
added the option to scale del2 and del4 mixing of u and theta based on mesh spacing<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/Registry
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2011-05-25 20:48:44 UTC (rev 852)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2011-05-25 21:04:15 UTC (rev 853)
@@ -6,6 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
+namelist logical sw_model config_h_ScaleWithMesh false
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_v_mom_eddy_visc2 0.0
@@ -72,6 +73,10 @@
var real zVertex ( nVertices ) iro zVertex - -
var integer indexToVertexID ( nVertices ) iro indexToVertexID - -
+var real meshDensity ( nCells ) iro meshDensity - -
+var real meshScalingDel2 ( nEdges ) ro meshScalingDel2 - -
+var real meshScalingDel4 ( nEdges ) ro meshScalingDel4 - -
+
var integer cellsOnEdge ( TWO nEdges ) iro cellsOnEdge - -
var integer nEdgesOnCell ( nCells ) iro nEdgesOnCell - -
var integer nEdgesOnEdge ( nEdges ) iro nEdgesOnEdge - -
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2011-05-25 20:48:44 UTC (rev 852)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2011-05-25 21:04:15 UTC (rev 853)
@@ -542,6 +542,7 @@
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, h_theta_eddy_visc4
real (kind=RKIND) :: u_diffusion
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, p_s
+ real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, divergence, vorticity, ke, pv_edge, geopotential, theta, ww, &
h_diabatic, tend_theta, theta_phys_tend, u_phys_tend
@@ -617,6 +618,10 @@
nVertLevels = grid % nVertLevels
nVertices = grid % nVertices
+
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
h_mom_eddy_visc2 = config_h_mom_eddy_visc2
h_mom_eddy_visc4 = config_h_mom_eddy_visc4
v_mom_eddy_visc2 = config_v_mom_eddy_visc2
@@ -699,7 +704,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 = meshScalingDel2(iEdge) * h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
end do
@@ -779,7 +784,8 @@
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
+ u_diffusion = meshScalingDel4(iEdge) * h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
end do
end do
@@ -867,7 +873,7 @@
cell2 = grid % cellsOnEdge % array(2,iEdge)
do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = meshScalingDel2(iEdge)*h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
tend_theta(k,cell1) = tend_theta(k,cell1) + flux
tend_theta(k,cell2) = tend_theta(k,cell2) - flux
@@ -906,7 +912,7 @@
cell2 = grid % cellsOnEdge % array(2,iEdge)
do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = meshScalingDel4(iEdge)*h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/mpas_interface.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/mpas_interface.F 2011-05-25 20:48:44 UTC (rev 852)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/mpas_interface.F 2011-05-25 21:04:15 UTC (rev 853)
@@ -29,6 +29,7 @@
call compute_solver_constants(block % time_levs(1) % state, mesh)
call compute_state_diagnostics(block % time_levs(1) % state, mesh)
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
+ call compute_mesh_scaling(mesh)
call initialize_advection_rk(mesh)
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
@@ -72,3 +73,37 @@
implicit none
end subroutine mpas_finalize
+
+
+subroutine compute_mesh_scaling(mesh)
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (grid_meta), 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
+
</font>
</pre>