<p><b>ringler@lanl.gov</b> 2011-03-01 10:19:11 -0700 (Tue, 01 Mar 2011)</p><p><br>
added high-order reconstruction for vertical transport.<br>
<br>
the high-order reconstruction is obtained by setting<br>
config_theta_vadv_order = {2,3,4}<br>
config_scalar_vadv_order = {2,3,4}<br>
<br>
when the value == 2, the new code reproduces bit-for-bit the previous version for the tests conducted.<br>
<br>
testing was done with test case # 2, baroclinic eddy with perturbation on a global 60 km mesh.<br>
<br>
the namelist.input.hyd_atmos is configured with values appropriate for a 60 km aqua-planet simulation.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/namelist.input.hyd_atmos
===================================================================
--- branches/mpas_cam_coupling/namelist.input.hyd_atmos 2011-02-28 18:45:22 UTC (rev 747)
+++ branches/mpas_cam_coupling/namelist.input.hyd_atmos 2011-03-01 17:19:11 UTC (rev 748)
@@ -1,20 +1,24 @@
&sw_model
config_test_case = 2
config_time_integration = 'SRK3'
- config_dt = 3600
- config_ntimesteps = 240
- config_output_interval = 24
+ config_dt = 225.0
+ config_ntimesteps = 11520
+ config_output_interval = 11520
config_number_of_sub_steps = 4
config_h_mom_eddy_visc2 = 0.0
- config_h_mom_eddy_visc4 = 0.0
+ config_h_mom_eddy_visc4 = 5.0e13
config_v_mom_eddy_visc2 = 0.0
config_h_theta_eddy_visc2 = 0.0
- config_h_theta_eddy_visc4 = 0.0
+ config_h_theta_eddy_visc4 = 5.0e13
config_v_theta_eddy_visc2 = 0.0
- config_theta_adv_order = 2
- config_scalar_adv_order = 2
+ config_theta_adv_order = 3
+ config_scalar_adv_order = 3
+ config_theta_vadv_order = 3
+ config_scalar_vadv_order = 3
+ config_apvm_upwinding = 0.5
+ config_positive_definite = .false.
+ config_monotonic = .true.
config_mp_physics = 0
- config_apvm_upwinding = 0.5
/
&io
Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/Registry
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2011-02-28 18:45:22 UTC (rev 747)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2011-03-01 17:19:11 UTC (rev 748)
@@ -13,12 +13,15 @@
namelist real sw_model config_h_theta_eddy_visc4 0.0
namelist real sw_model config_v_theta_eddy_visc2 0.0
namelist integer sw_model config_number_of_sub_steps 4
-namelist integer sw_model config_theta_adv_order 2
-namelist integer sw_model config_scalar_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic true
namelist integer sw_model config_mp_physics 0
namelist real sw_model config_apvm_upwinding 0.5
+namelist real sw_model config_coef_3rd_order 1.0
+namelist integer sw_model config_theta_vadv_order 2
+namelist integer sw_model config_scalar_vadv_order 2
+namelist integer sw_model config_theta_adv_order 2
+namelist integer sw_model config_scalar_adv_order 2
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/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2011-02-28 18:45:22 UTC (rev 747)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/module_time_integration.F 2011-03-01 17:19:11 UTC (rev 748)
@@ -559,6 +559,23 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND) :: coef_3rd_order
+
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+!-----------
+
+ coef_3rd_order = config_coef_3rd_order
+
+
h => s % h % array
u => s % u % array
h_edge => s % h_edge % array
@@ -988,11 +1005,37 @@
! Note: we are also dividing through by the cell area after the horizontal flux divergence
!
do iCell = 1, nCells
+
wdtn(1) = 0.
- do k=2,nVertLevels
+ if (config_theta_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
+ end do
+
+ else if (config_theta_vadv_order == 3) then
+
+ k = 2
wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
- end do
+ do k=3,nVertLevels-1
+ wdtn(k) = flux3( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), ww(k,iCell), coef_3rd_order )
+ end do
+ k = nVertLevels
+ wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
+
+ else if (config_theta_vadv_order == 4) then
+
+ k = 2
+ wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
+ do k=3,nVertLevels-1
+ wdtn(k) = flux4( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), ww(k,iCell) )
+ end do
+ k = nVertLevels
+ wdtn(k) = ww(k,icell)*(fnm(k)*theta(k,iCell)+fnp(k)*theta(k-1,iCell))
+
+ end if
wdtn(nVertLevels+1) = 0.
+
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdnw(k)*(wdtn(k+1)-wdtn(k))
!! tend_theta(k,iCell) = tend_theta(k) + h(k,iCell)*h_diabatic(k,iCell)
@@ -1325,6 +1368,16 @@
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
coef_3rd_order = 0.
if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
@@ -1460,14 +1513,41 @@
!
do iCell=1,grid % nCells
+ wdtn(:,1) = 0.
+ if (config_scalar_vadv_order == 2) then
+ do k = 2, nVertLevels
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ end do
+ end do
+ else
+ k = 2
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+ if ( config_scalar_vadv_order == 3 ) then
+ do k=3,nVertLevels-1
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = flux3( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), &
+ wwAvg(k,iCell), coef_3rd_order )
+ end do
+ end do
+ else if ( config_scalar_vadv_order == 4 ) then
+ do k=3,nVertLevels-1
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = flux4( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), wwAvg(k,iCell) )
+ end do
+ end do
+ end if
+ k = nVertLevels
+ do iScalar=1,num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
- wdtn(:,1) = 0.
- do k = 2, nVertLevels
- do iScalar=1,num_scalars
- wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- end do
- end do
- wdtn(:,nVertLevels+1) = 0.
+ end if
+ wdtn(:,nVertLevels+1) = 0.
do k=1,grid % nVertLevelsSolve
do iScalar=1,num_scalars
@@ -1533,6 +1613,16 @@
real (kind=RKIND), parameter :: eps=1.e-20
real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
+ ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+
+ flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
+ flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
scalar_phys_tend => s_old % scalars_phys_tend % array
@@ -1579,11 +1669,30 @@
do iCell=1,grid % nCells
if (k < grid % nVertLevels) then
+
+ if ((config_scalar_vadv_order == 2) .or. (k==1) .or. (k==grid % nVertLevels-1)) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * &
+ (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
+ end do
+ else if (config_scalar_vadv_order == 3) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * flux3( scalar_new(iScalar,k-1,iCell),scalar_new(iScalar,k ,iCell), &
+ scalar_new(iScalar,k+1,iCell),scalar_new(iScalar,k+2,iCell), &
+ wwAvg(k+1,iCell), coef_3rd_order )
+ end do
+ else if (config_scalar_vadv_order == 4) then
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * flux4( scalar_new(iScalar,k-1,iCell),scalar_new(iScalar,k ,iCell), &
+ scalar_new(iScalar,k+1,iCell),scalar_new(iScalar,k+2,iCell), wwAvg(k+1,iCell) )
+ end do
+ end if
+
+
cell_upwind = k
if (wwAvg(k+1,iCell) >= 0) cell_upwind = k+1
+
do iScalar=1,num_scalars
- v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * &
- (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
v_flux_upwind(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * scalar_old(iScalar,cell_upwind,iCell)
v_flux(iScalar,iCell,km0) = v_flux(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km0)
! v_flux(iScalar,iCell,km0) = 0. ! use only upwind - for testing
</font>
</pre>