<p><b>duda</b> 2012-05-18 14:56:37 -0600 (Fri, 18 May 2012)</p><p>BRANCH COMMIT<br>
<br>
- Correct logic for applying "2d_fixed" or "2d_smagorinsky" horizontal mixing<br>
<br>
- Add trap to avoid small edge lengths (<0.25*dcEdge) in horizontal momentum mixing<br>
<br>
- Corrections to use of Prandtl number (K_m = prandtl_inv * K_v (prandtl_inv, computed from prandtl specified in framework/mpas_constants.F))<br>
<br>
- Compute derived constants in mpas_constants.F, rather than specifying numerical values<br>
<br>
<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
M src/framework/mpas_constants.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-05-18 18:54:23 UTC (rev 1920)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-05-18 20:56:37 UTC (rev 1921)
@@ -169,7 +169,7 @@
if (debug) write(0,*) ' compute_dyn_tend '
block => domain % blocklist
do while (associated(block))
- call atm_compute_dyn_tend( block % tend, block % state % time_levs(2) % state, block % diag, block % mesh, rk_step )
+ call atm_compute_dyn_tend( block % tend, block % state % time_levs(2) % state, block % diag, block % mesh, rk_step, dt )
block => block % next
end do
if (debug) write(0,*) ' finished compute_dyn_tend '
@@ -1792,7 +1792,7 @@
!----
- subroutine atm_compute_dyn_tend(tend, s, diag, grid, rk_step)
+ subroutine atm_compute_dyn_tend(tend, s, diag, grid, rk_step, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -1811,6 +1811,7 @@
type (diag_type), intent(in) :: diag
type (mesh_type), intent(in) :: grid
integer, intent(in) :: rk_step
+ real (kind=RKIND), intent(in) :: dt
logical, parameter :: rk_diffusion = .false.
@@ -1854,6 +1855,7 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
real (kind=RKIND) :: cf1, cf2, cf3, pr, pl
+ real (kind=RKIND) :: prandtl_inv
! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
@@ -1981,6 +1983,8 @@
adv_coefs => grid % adv_coefs % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
+ prandtl_inv = 1.0_RKIND/prandtl
+
!
! Compute u (normal) velocity tendency for each edge (cell face)
!
@@ -2003,6 +2007,7 @@
end do
do k=1, nVertLevels
kdiff(k,iCell) = (c_s * config_len_disp)**2 * sqrt(d_diag(k)**2 + d_off_diag(k)**2)
+ kdiff(k,iCell) = min(kdiff(k,iCell),(0.01*config_len_disp**2)/dt)
end do
end do
delsq_horiz_mixing = .true.
@@ -2166,7 +2171,7 @@
if (delsq_horiz_mixing) then
- if (h_mom_eddy_visc2 > 0.0) then
+ if ((h_mom_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2180,7 +2185,7 @@
! only valid for h_mom_eddy_visc2 == 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))
u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
u_diffusion = u_diffusion * meshScalingDel2(iEdge)
@@ -2203,7 +2208,7 @@
! only valid for h_mom_eddy_visc2 == 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))
u_diffusion = rho_edge(k,iEdge)* 0.5*(kdiff(k,cell1)+kdiff(k,cell2)) * u_diffusion
u_diffusion = u_diffusion * meshScalingDel2(iEdge)
@@ -2215,7 +2220,7 @@
end if ! delsq_horiz_mixing for u
- if ( h_mom_eddy_visc4 > 0.0 ) then
+ if ((h_mom_eddy_visc4 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_u(nVertLevels, nEdges+1))
@@ -2525,7 +2530,7 @@
if (delsq_horiz_mixing) then
- if (h_mom_eddy_visc2 > 0.0) then
+ if ((h_mom_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
@@ -2569,7 +2574,7 @@
end if
- if ( h_mom_eddy_visc4 > 0.0 ) then
+ if ( (h_mom_eddy_visc4 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
allocate(delsq_theta(nVertLevels, nCells+1))
@@ -2820,7 +2825,7 @@
tend_theta_euler = 0.
if (delsq_horiz_mixing) then
- if ( h_theta_eddy_visc2 > 0.0 ) then
+ if ( (h_theta_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed") ) then
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
@@ -2828,7 +2833,7 @@
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
+ theta_turb_flux = h_theta_eddy_visc2*prandtl_inv*(theta_m(k,cell2) - theta_m(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
@@ -2848,7 +2853,7 @@
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,grid % nVertLevels
- theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl &
+ theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl_inv &
*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
@@ -2864,7 +2869,7 @@
end if
- if ( h_theta_eddy_visc4 > 0.0 ) then
+ if ( (h_theta_eddy_visc4 > 0.0) .and. (config_horiz_mixing == "2d_fixed") ) then
allocate(delsq_theta(nVertLevels, nCells+1))
@@ -2892,7 +2897,7 @@
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
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 = h_theta_eddy_visc4*prandtl_inv*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel4(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
@@ -2976,10 +2981,10 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
-! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl_inv*rho_zz(k,iCell)*(&
! (theta_m(k+1,iCell)-theta_m(k ,iCell))/(zp-z0) &
! -(theta_m(k ,iCell)-theta_m(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
- tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
+ tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl_inv*rho_zz(k,iCell)*(&
(theta_m(k+1,iCell)-theta_m(k ,iCell))/(zp-z0) &
-(theta_m(k ,iCell)-theta_m(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -2998,10 +3003,10 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
-! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl_inv*rho_zz(k,iCell)*(&
! ((theta_m(k+1,iCell)-t_init(k+1,iCell))-(theta_m(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
! -((theta_m(k ,iCell)-t_init(k,iCell))-(theta_m(k-1,iCell)-t_init(k-1,iCell)))/(z0-zm) )/(0.5*(zp-zm))
- tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
+ tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl_inv*rho_zz(k,iCell)*(&
((theta_m(k+1,iCell)-t_init(k+1,iCell))-(theta_m(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
-((theta_m(k ,iCell)-t_init(k,iCell))-(theta_m(k-1,iCell)-t_init(k-1,iCell)))/(z0-zm) )/(0.5*(zp-zm))
end do
Modified: branches/atmos_physics/src/framework/mpas_constants.F
===================================================================
--- branches/atmos_physics/src/framework/mpas_constants.F 2012-05-18 18:54:23 UTC (rev 1920)
+++ branches/atmos_physics/src/framework/mpas_constants.F 2012-05-18 20:56:37 UTC (rev 1921)
@@ -8,8 +8,8 @@
real (kind=RKIND), parameter :: gravity = 9.80616
real (kind=RKIND), parameter :: rgas = 287.
real (kind=RKIND), parameter :: cp = 1003.
- real (kind=RKIND), parameter :: cv = 716. ! cp - rgas
- real (kind=RKIND), parameter :: cvpm = -.71385842 ! -cv/cp
+ real (kind=RKIND), parameter :: cv = cp - rgas
+ real (kind=RKIND), parameter :: cvpm = -cv/cp
real (kind=RKIND), parameter :: prandtl = 1.0
</font>
</pre>