<p><b>duda</b> 2011-06-10 14:06:17 -0600 (Fri, 10 Jun 2011)</p><p>BRANCH COMMIT<br>
<br>
Add code to scale horizontal diffusion coefficients with mesh resolution, <br>
as is done in the SW model of the repository trunk. There are currently <br>
two differences from the model trunk:<br>
<br>
1) The mesh density field is called densityFunction in the input and output<br>
files, but is internally named meshDensity, since we have a set of existing<br>
mesh files that we have yet to update (by changing densityFunction to<br>
meshDensity)<br>
<br>
2) The del2 viscosities are scaled with to the 0.5 power of the mesh density<br>
(rather than the 5/12 power), and the del4 viscosities are scaled directly<br>
by the mesh density (rather than the 5/6 power).<br>
<br>
<br>
Also, write out min/max of u, w, and scalars as the global min/max to aid in <br>
testing; local min/max values can be written out by uncommenting a few lines<br>
in module_time_integration.F<br>
<br>
<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_mpas_core.F<br>
M src/core_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_time_integration.F<br>
M namelist.input.nhyd_atmos_jw<br>
M namelist.input.nhyd_atmos<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/namelist.input.nhyd_atmos
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/namelist.input.nhyd_atmos 2011-06-10 20:06:17 UTC (rev 898)
@@ -17,6 +17,7 @@
config_monotonic = .false.
config_epssm = 0.1
config_smdiv = 0.1
+ config_h_ScaleWithMesh = .false.
/
&dimensions
Modified: branches/atmos_physics/namelist.input.nhyd_atmos_jw
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos_jw 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/namelist.input.nhyd_atmos_jw 2011-06-10 20:06:17 UTC (rev 898)
@@ -25,6 +25,7 @@
config_monotonic = .false.
config_epssm = 0.1
config_smdiv = 0.1
+ config_h_ScaleWithMesh = .false.
/
&dimensions
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-06-10 20:06:17 UTC (rev 898)
@@ -97,6 +97,8 @@
var persistent real kiteAreasOnVertex ( vertexDegree nVertices ) 0 io kiteAreasOnVertex mesh - -
var persistent real fEdge ( nEdges ) 0 io fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 io fVertex mesh - -
+
+var persistent real densityFunction ( nCells ) 0 iro densityFunction mesh - -
# some solver scalar coefficients
Modified: branches/atmos_physics/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-06-10 20:06:17 UTC (rev 898)
@@ -31,6 +31,7 @@
namelist real nhyd_model config_epssm 0.1
namelist real nhyd_model config_smdiv 0.1
namelist real nhyd_model config_apvm_upwinding 0.5
+namelist logical nhyd_model config_h_ScaleWithMesh false
namelist integer dimensions config_nvertlevels 26
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -111,6 +112,10 @@
var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
+var persistent real densityFunction ( nCells ) 0 iro meshDensity mesh - -
+var persistent real meshScalingDel2 ( nEdges ) 0 ro meshScalingDel2 mesh - -
+var persistent real meshScalingDel4 ( nEdges ) 0 ro meshScalingDel4 mesh - -
+
# some solver scalar coefficients
# coefficients for vertical extrapolation to the surface
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-06-10 20:06:17 UTC (rev 898)
@@ -111,6 +111,13 @@
#endif
if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
+
+ call compute_mesh_scaling(mesh)
+
+ write(0,*) 'min/max of meshScalingDel2 = ', minval(mesh % meshScalingDel2 % array(1:mesh%nEdges)), &
+ maxval(mesh % meshScalingDel2 % array(1:mesh%nEdges))
+ write(0,*) 'min/max of meshScalingDel4 = ', minval(mesh % meshScalingDel4 % array(1:mesh%nEdges)), &
+ maxval(mesh % meshScalingDel4 % array(1:mesh%nEdges))
end subroutine mpas_init_block
@@ -264,4 +271,36 @@
end subroutine mpas_core_finalize
+
+ subroutine compute_mesh_scaling(mesh)
+
+ use grid_types
+
+ 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)**0.5
+ meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)
+ end do
+ end if
+
+ end subroutine compute_mesh_scaling
+
end module mpas_core
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-06-10 19:03:45 UTC (rev 897)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-06-10 20:06:17 UTC (rev 898)
@@ -410,7 +410,8 @@
! end if
! if(debug) then
- 101 format(' min, max scalar',i4,2(1x,e17.10))
+ 101 format('local min, max scalar',i4,2(1x,e17.10))
+ 102 format('global min, max scalar',i4,2(1x,e17.10))
write(0,*)
block => domain % blocklist
do while (associated(block))
@@ -422,7 +423,10 @@
scalar_max = max(scalar_max, block % state % time_levs(2) % state % w % array(k,iCell))
enddo
enddo
- write(0,*) ' min, max w ',scalar_min, scalar_max
+ call dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
+ call dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
+! write(0,*) 'local min, max w ',scalar_min, scalar_max
+ write(0,*) 'global min, max w ',global_scalar_min, global_scalar_max
scalar_min = 0.
scalar_max = 0.
@@ -432,7 +436,10 @@
scalar_max = max(scalar_max, block % state % time_levs(2) % state % u % array(k,iEdge))
enddo
enddo
- write(0,*) ' min, max u ',scalar_min, scalar_max
+ call dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
+ call dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
+! write(0,*) 'local min, max u ',scalar_min, scalar_max
+ write(0,*) 'global min, max u ',global_scalar_min, global_scalar_max
do iScalar = 1, block % state % time_levs(2) % state % num_scalars
scalar_min = 0.
@@ -443,7 +450,10 @@
scalar_max = max(scalar_max, block % state % time_levs(2) % state % scalars % array(iScalar,k,iCell))
enddo
enddo
- write(0,101) iScalar,scalar_min,scalar_max
+ call dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
+ call dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
+! write(0,101) iScalar,scalar_min,scalar_max
+ write(0,102) iScalar,global_scalar_min,global_scalar_max
end do
block => block % next
@@ -1182,7 +1192,7 @@
real (kind=RKIND), dimension( s_old % num_scalars, grid % nVertLevels + 1 ) :: wdtn
integer :: nVertLevels
- real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
+ real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, scalar_turb_flux, z1,z2,z3,z4,zm,z0,zp
@@ -1226,6 +1236,8 @@
fnm => grid % fzm % array
fnp => grid % fzp % array
rdnw => grid % rdzw % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
nVertLevels = grid % nVertLevels
@@ -1358,6 +1370,7 @@
do iScalar=1,s_old % num_scalars
scalar_turb_flux = h_theta_eddy_visc2*prandtl* &
(scalar_new(iScalar,k,cell2) - scalar_new(iScalar,k,cell1))/dcEdge(iEdge)
+ scalar_turb_flux = scalar_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * scalar_turb_flux
scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) + flux/areaCell(cell1)
scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) - flux/areaCell(cell2)
@@ -1378,6 +1391,7 @@
do iScalar=1,s_old % num_scalars
scalar_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl* &
(scalar_new(iScalar,k,cell2) - scalar_new(iScalar,k,cell1))/dcEdge(iEdge)
+ scalar_turb_flux = scalar_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * scalar_turb_flux
scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) + flux/areaCell(cell1)
scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) - flux/areaCell(cell2)
@@ -1922,7 +1936,7 @@
real (kind=RKIND) :: h_mom_eddy_visc2, v_mom_eddy_visc2, h_mom_eddy_visc4
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 :: fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:), pointer :: fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, zgrid, rho_edge, rho, ru, u, v, tend_u, &
circulation, divergence, vorticity, ke, pv_edge, theta, rw, tend_rho, &
rt_diabatic_tend, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
@@ -2012,7 +2026,10 @@
defc_a => grid % defc_a % array
defc_b => grid % defc_b % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
tend_u => tend % u % array
tend_theta => tend % theta % array
tend_w => tend % w % array
@@ -2216,6 +2233,7 @@
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
+ u_diffusion = u_diffusion * meshScalingDel2(iEdge)
! tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + u_diffusion
@@ -2238,6 +2256,7 @@
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = rho_edge(k,iEdge)* 0.5*(kdiff(k,cell1)+kdiff(k,cell2)) * u_diffusion
+ u_diffusion = u_diffusion * meshScalingDel2(iEdge)
! tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + u_diffusion
@@ -2269,7 +2288,7 @@
! only valid for h_mom_eddy_visc4 == constant
!
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / max(dvEdge(iEdge), 0.25*dcEdge(iEdge))
delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
end do
@@ -2318,10 +2337,11 @@
! only valid for h_mom_eddy_visc4 == constant
!
u_diffusion = rho_edge(k,iEdge) * ( ( delsq_divergence(k,cell2) - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -( delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge) &
+ -( delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) / max(dvEdge(iEdge), 0.25*dcEdge(iEdge)) &
)
! tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ u_diffusion = u_diffusion * meshScalingDel4(iEdge)
tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
end do
end do
@@ -2527,6 +2547,7 @@
do k=2,grid % nVertLevels
theta_turb_flux = h_mom_eddy_visc2*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
! tend_w(k,cell1) = tend_w(k,cell1) + flux
! tend_w(k,cell2) = tend_w(k,cell2) - flux
@@ -2547,6 +2568,7 @@
! theta_turb_flux = h_mom_eddy_visc2*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
theta_turb_flux = 0.25*(kdiff(k,cell1)+kdiff(k,cell2)+kdiff(k-1,cell1)+kdiff(k,cell2)) &
*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
! tend_w(k,cell1) = tend_w(k,cell1) + flux
! tend_w(k,cell2) = tend_w(k,cell2) - flux
@@ -2593,6 +2615,7 @@
do k=2,grid % nVertLevels
theta_turb_flux = h_mom_eddy_visc4*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = theta_turb_flux * meshScalingDel4(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
! tend_w(k,cell1) = tend_w(k,cell1) - flux
@@ -2838,6 +2861,7 @@
do k=1,grid % nVertLevels
theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
! tend_theta(k,cell1) = tend_theta(k,cell1) + flux
! tend_theta(k,cell2) = tend_theta(k,cell2) - flux
@@ -2858,6 +2882,7 @@
do k=1,grid % nVertLevels
theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl &
*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
! tend_theta(k,cell1) = tend_theta(k,cell1) + flux
! tend_theta(k,cell2) = tend_theta(k,cell2) - flux
@@ -2900,6 +2925,7 @@
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 = theta_turb_flux * meshScalingDel4(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
! tend_theta(k,cell1) = tend_theta(k,cell1) - flux
</font>
</pre>