<p><b>duda</b> 2010-07-28 15:30:42 -0600 (Wed, 28 Jul 2010)</p><p>BRANCH COMMIT<br>
<br>
New additions to the nonhydrostatic core from Bill:<br>
- 2d Smagorinsky horizontal mixing<br>
- higher-order vertical advection for scalars<br>
- a new namelist for running the J&W baroclinic wave case<br>
using the non-hydrostatic core<br>
- accompanying changes to Registry and namelists<br>
<br>
<br>
M namelist.input.nhyd_atmos_squall<br>
M graphics/ncl/cells_nhyd_sph1.ncl<br>
M src/core_nhyd_atmos/mpas_interface.F<br>
M src/core_nhyd_atmos/module_advection.F<br>
M src/core_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_time_integration.F<br>
M namelist.input.nhyd_atmos<br>
A namelist.input.nhyd_atmos_jw<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/graphics/ncl/cells_nhyd_sph1.ncl
===================================================================
--- branches/atmos_nonhydrostatic/graphics/ncl/cells_nhyd_sph1.ncl 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/graphics/ncl/cells_nhyd_sph1.ncl 2010-07-28 21:30:42 UTC (rev 444)
@@ -90,7 +90,12 @@
res@cnLineLabelsOn = True
end if
-; res@cnLevelSpacingF = 10.0
+; res@cnLevelSelectionMode = 2
+; res@cnLevels = (/950.,960.,970.,980.,990.,1000.,1010.,1020./)
+ res@cnLevelSelectionMode = 3
+ res@cnLevelSpacingF = 2.
+ res@cnMinLevelValF = 940.
+ res@cnMaxLevelValF = 1020.
res@cnInfoLabelOn = True
res@lbLabelAutoStride = True
@@ -108,8 +113,6 @@
res@mpPerimOn = False
res@gsnFrame = False
- res@cnLevelSelectionMode = 2
- res@cnLevels = (/950.,960.,970.,980.,990.,1000.,1010.,1020./)
t = stringtointeger(getenv("T"))
if (plotfield .eq. "h") then
@@ -155,6 +158,8 @@
; fld = f->pressure(t,:,25)+f->pressure_base(:,25)
+; fld = f->kdiff(t,:,0)
+
; zg = f->zgrid
; csizes = dimsizes(pfirst)
; fld = pfirst
Modified: branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos 2010-07-28 21:30:42 UTC (rev 444)
@@ -1,4 +1,4 @@
-&sw_model
+&nhyd_model
config_test_case = 2
config_time_integration = 'SRK3'
config_dt = 1800
Added: branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_jw
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_jw (rev 0)
+++ branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_jw 2010-07-28 21:30:42 UTC (rev 444)
@@ -0,0 +1,40 @@
+&nhyd_model
+ config_test_case = 2
+ config_time_integration = 'SRK3'
+ config_dt = 450
+ config_ntimesteps = 1920
+ config_output_interval = 192
+ config_number_of_sub_steps = 6
+ config_h_mom_eddy_visc2 = 0.0e+04
+ config_h_mom_eddy_visc4 = 0.
+ config_v_mom_eddy_visc2 = 00.0
+ config_h_theta_eddy_visc2 = 0.0e+04
+ config_h_theta_eddy_visc4 = 00.
+ config_v_theta_eddy_visc2 = 00.0
+ config_horiz_mixing = '2d_smagorinsky'
+ config_len_disp = 60000.
+ config_u_vadv_order = 3
+ config_w_vadv_order = 3
+ config_theta_vadv_order = 3
+ config_scalar_vadv_order = 3
+ config_theta_adv_order = 3
+ config_scalar_adv_order = 3
+ config_scalar_advection = .false.
+ config_positive_definite = .false.
+ config_coef_3rd_order = 1.0
+ config_monotonic = .false.
+ config_epssm = 0.1
+ config_smdiv = 0.1
+/
+
+&io
+ config_input_name = 'grid.nc'
+ config_output_name = 'output.nc'
+ config_restart_name = 'restart.nc'
+/
+
+&restart
+ config_restart_interval = 3000
+ config_do_restart = .false.
+ config_restart_time = 1036800.0
+/
Modified: branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_squall
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_squall 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_squall 2010-07-28 21:30:42 UTC (rev 444)
@@ -1,4 +1,4 @@
-&sw_model
+&nhyd_model
config_test_case = 4
config_time_integration = 'SRK3'
config_dt = 6.
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-07-28 21:30:42 UTC (rev 444)
@@ -1,25 +1,35 @@
#
# namelist type namelist_record name default_value
#
-namelist integer sw_model config_test_case 5
-namelist character sw_model config_time_integration SRK3
-namelist real sw_model config_dt 172.8
-namelist integer sw_model config_ntimesteps 7500
-namelist integer sw_model config_output_interval 500
-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
-namelist real sw_model config_h_theta_eddy_visc2 0.0
-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_epssm 0.1
-namelist real sw_model config_smdiv 0.1
+namelist integer nhyd_model config_test_case 5
+namelist character nhyd_model config_time_integration SRK3
+namelist real nhyd_model config_dt 172.8
+namelist integer nhyd_model config_ntimesteps 7500
+namelist integer nhyd_model config_output_interval 500
+namelist character nhyd_model config_horiz_mixing 2d_smagorinsky
+namelist real nhyd_model config_h_mom_eddy_visc2 0.0
+namelist real nhyd_model config_h_mom_eddy_visc4 0.0
+namelist real nhyd_model config_v_mom_eddy_visc2 0.0
+namelist real nhyd_model config_h_theta_eddy_visc2 0.0
+namelist real nhyd_model config_h_theta_eddy_visc4 0.0
+namelist real nhyd_model config_v_theta_eddy_visc2 0.0
+namelist integer nhyd_model config_number_of_sub_steps 4
+namelist integer nhyd_model config_w_adv_order 2
+namelist integer nhyd_model config_theta_adv_order 2
+namelist integer nhyd_model config_scalar_adv_order 2
+namelist integer nhyd_model config_u_vadv_order 2
+namelist integer nhyd_model config_w_vadv_order 2
+namelist integer nhyd_model config_theta_vadv_order 2
+namelist integer nhyd_model config_scalar_vadv_order 2
+namelist real nhyd_model config_coef_3rd_order 1.0
+namelist logical nhyd_model config_scalar_advection true
+namelist logical nhyd_model config_positive_definite false
+namelist logical nhyd_model config_monotonic true
+namelist logical nhyd_model config_mix_full true
+namelist real nhyd_model config_len_disp 0.
+namelist integer nhyd_model config_mp_physics 0.
+namelist real nhyd_model config_epssm 0.1
+namelist real nhyd_model config_smdiv 0.1
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -220,6 +230,11 @@
var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
var integer advCells ( TWENTYONE nCells ) - advCells - -
+# Space needed for deformation calculation weights
+var real defc_a ( maxEdges nCells ) - defc_a - -
+var real defc_b ( maxEdges nCells ) - defc_b - -
+var real kdiff ( nVertLevels nCells Time ) - kdiff - -
+
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F 2010-07-28 21:30:42 UTC (rev 444)
@@ -685,4 +685,194 @@
!
END SUBROUTINE ELGS
+!-------------------------------------------------------------
+
+ subroutine initialize_deformation_weights( grid )
+
+!
+! compute the cell coefficients for the deformation calculations
+! WCS, 13 July 2010
+!
+ implicit none
+
+ type (grid_meta), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+
+! local variables
+
+ real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ real (kind=RKIND), dimension(grid % nEdges) :: xe, ye
+ real (kind=RKIND), dimension(grid % nCells) :: theta_abs
+
+ real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
+ real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
+ real (kind=RKIND) :: xm, ym, zm, dl, xec, yec, zec
+ real (kind=RKIND) :: thetae_tmp, xe_tmp, ye_tmp
+ real (kind=RKIND) :: xv1, xv2, yv1, yv2, zv1, zv2
+ integer :: i, j, k, ip1, ip2, m, n, ip1a, ii
+ integer :: iCell, iEdge
+ real (kind=RKIND) :: pii
+ real (kind=RKIND) :: x0, y0, x1, y1, x2, y2, x3, y3, x4, y4, x5, y5
+ real (kind=RKIND) :: pdx1, pdx2, pdx3, pdy1, pdy2, pdy3, dx1, dx2, dy1, dy2
+ real (kind=RKIND) :: angv1, angv2, dl1, dl2
+ real (kind=RKIND), dimension(25) :: dxe, dye, x2v, y2v, xp, yp, xpt, ypt
+
+ real (kind=RKIND) :: length_scale
+ integer :: ma,na, cell_add, mw, nn
+ integer, dimension(25) :: cell_list
+
+ integer :: cell1, cell2, iv
+ logical :: do_the_cell
+ real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, sumw1, sumw2, xptt, area_cellt
+
+ logical, parameter :: debug = .false.
+
+ if (debug) write(0,*) ' in def weight calc '
+
+ defc_a => grid % defc_a % array
+ defc_b => grid % defc_b % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+
+ defc_a(:,:) = 0.
+ defc_b(:,:) = 0.
+
+ pii = 2.*asin(1.0)
+
+ if (debug) write(0,*) ' beginning cell loop '
+
+ do iCell = 1, grid % nCells
+
+ if (debug) write(0,*) ' cell loop ', iCell
+
+ cell_list(1) = iCell
+ do i=2, grid % nEdgesOnCell % array(iCell)+1
+ cell_list(i) = grid % CellsOnCell % array(i-1,iCell)
+ end do
+ n = grid % nEdgesOnCell % array(iCell) + 1
+
+! check to see if we are reaching outside the halo
+
+ if (debug) write(0,*) ' points ', n
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
+ end do
+
+
+ if ( .not. do_the_cell ) cycle
+
+
+! compute poynomial fit for this cell if all needed neighbors exist
+ if ( grid % on_a_sphere ) then
+
+ xc(1) = grid % xCell % array(iCell)/a
+ yc(1) = grid % yCell % array(iCell)/a
+ zc(1) = grid % zCell % array(iCell)/a
+
+
+ do i=2,n
+ iv = grid % verticesOnCell % array(i-1,iCell)
+ xc(i) = grid % xVertex % array(iv)/a
+ yc(i) = grid % yVertex % array(iv)/a
+ zc(i) = grid % zVertex % array(iv)/a
+ end do
+
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
+
+! angles from cell center to neighbor centers (thetav)
+
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) ip2 = 2
+
+ thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1), &
+ xc(ip2), yc(ip2), zc(ip2) )
+
+ dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
+ end do
+
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ end do
+
+ thetat(1) = 0. ! this defines the x direction, cell center 1 ->
+! thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
+ do i=2,n-1
+ thetat(i) = thetat(i-1) + thetav(i-1)
+ end do
+
+ do i=1,n-1
+ xp(i) = cos(thetat(i)) * dl_sphere(i)
+ yp(i) = sin(thetat(i)) * dl_sphere(i)
+ end do
+
+ else ! On an x-y plane
+
+ xp(1) = grid % xCell % array(iCell)
+ yp(1) = grid % yCell % array(iCell)
+
+
+ do i=2,n
+ iv = grid % verticesOnCell % array(i-1,iCell)
+ xp(i) = grid % xVertex % array(iv)
+ yp(i) = grid % yVertex % array(iv)
+ end do
+
+ end if
+
+! thetat(1) = 0.
+ thetat(1) = theta_abs(iCell)
+ do i=2,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ thetat(i) = plane_angle( 0.,0.,0., &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
+ 0., 0., 1.)
+ thetat(i) = thetat(i) + thetat(i-1)
+ end do
+
+ area_cell = 0.
+ area_cellt = 0.
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
+ area_cellt = area_cellt + (0.25*(xp(i)+xp(ip1))*cos(thetat(i)) + 0.25*(yp(i)+yp(ip1))*sin(thetat(i)))*dl
+ end do
+ if (debug) write(0,*) ' area_cell, area_cellt ',area_cell, area_cellt,area_cell-area_cellt
+
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ sint2 = (sin(thetat(i)))**2
+ cost2 = (cos(thetat(i)))**2
+ sint_cost = sin(thetat(i))*cos(thetat(i))
+ defc_a(i,iCell) = dl*(cost2 - sint2)/area_cell
+ defc_b(i,iCell) = dl*2.*sint_cost/area_cell
+ if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
+ defc_a(i,iCell) = - defc_a(i,iCell)
+ defc_b(i,iCell) = - defc_b(i,iCell)
+ end if
+
+ end do
+
+ end do
+
+ if (debug) write(0,*) ' exiting def weight calc '
+
+ end subroutine initialize_deformation_weights
+
end module advection
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-07-28 21:30:42 UTC (rev 444)
@@ -73,8 +73,8 @@
logical, parameter :: debug = .false.
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
- logical, parameter :: do_microphysics = .true.
- logical, parameter :: scalar_advection = .true.
+! logical, parameter :: do_microphysics = .true.
+ logical, parameter :: do_microphysics = .false.
real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
real (kind=RKIND) :: global_domain_mass, global_scalar_mass, global_scalar_min, global_scalar_max
@@ -248,7 +248,7 @@
! ************ advection of moist variables here...
- if(scalar_advection) then
+ if (config_scalar_advection) then
block => domain % blocklist
do while (associated(block))
@@ -776,7 +776,7 @@
- rdzw(k)*(dpzx(k+1)-dpzx(k))
pgrad = 0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad
du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
-! + (0.005/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
+! + (0.05/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
ru_p(k,iEdge) = ru_p(k,iEdge) + du(k)
@@ -1077,7 +1077,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho, zgrid
+ real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho, zgrid, kdiff
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, qv_init
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -1087,17 +1087,27 @@
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
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
- logical, parameter :: mix_full = .false.
-! logical, parameter :: mix_full = .true.
- 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
+ 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*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
+
+ coef_3rd_order = config_coef_3rd_order
+
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
+ kdiff => s_new % kdiff % array
deriv_two => grid % deriv_two % array
!**** uhAvg => grid % uhAvg % array
uhAvg => grid % ruAvg % array
@@ -1188,17 +1198,6 @@
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
-! old version of the above code, with coef_3rd_order assumed to be 1.0
-! if (uhAvg(k,iEdge) > 0) then
-! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
-! else
-! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
-! end if
-
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)
@@ -1265,6 +1264,26 @@
end if
end do
+ else if ( config_horiz_mixing == "2d_smagorinsky") then
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+ do iScalar=1,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)
+ 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)
+ end do
+ end do
+
+ end if
+ end do
+
end if
! vertical mixing
@@ -1289,7 +1308,7 @@
end do
end do
- if ( .not. mix_full) then
+ if ( .not. config_mix_full) then
iScalar = index_qv
do k=2,nVertLevels-1
z1 = zgrid(k-1,iCell)
@@ -1317,14 +1336,42 @@
do iCell=1,grid % nCells
- 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.
+ 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
+ end if
+ wdtn(:,nVertLevels+1) = 0.
+
do k=1,grid % nVertLevelsSolve
do iScalar=1,num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
@@ -1378,6 +1425,18 @@
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*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+!------
+
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
deriv_two => grid % deriv_two % array
@@ -1416,10 +1475,12 @@
scale_out(:,:,:) = 1.
scale_in(:,:,:) = 1.
- 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
+! 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
+ coef_3rd_order = config_coef_3rd_order
+
do k = 1, grid % nVertLevels
kcp1 = min(k+1,grid % nVertLevels)
kcm1 = max(k-1,1)
@@ -1429,17 +1490,37 @@
do iCell=1,grid % nCells
if (k < grid % nVertLevels) then
+
+ if ((config_scalar_vadv_order == 2) .or. (k==1) .or. (k==grid % nVertLevels)) 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
s_update(iScalar,iCell,km0) = scalar_old(iScalar,k,iCell) * h_old(k,iCell) &
- rdnw(k) * (v_flux_upwind(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km1))
end do
+
+
else
do iScalar=1,num_scalars
v_flux(iScalar,iCell,km0) = 0.
@@ -1716,7 +1797,7 @@
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, &
h_diabatic, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
- h_divergence
+ h_divergence, kdiff
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -1736,17 +1817,36 @@
! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
- logical, parameter :: mix_full = .false.
-! logical, parameter :: mix_full = .true.
- integer :: w_adv_order
+ real (kind=RKIND), parameter :: c_s = 0.25
+ real (kind=RKIND), dimension( grid % nVertLevels ) :: d_diag, d_off_diag
+ real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
+ logical :: delsq_horiz_mixing
+
+
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*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+
+!-----------
+
+ coef_3rd_order = config_coef_3rd_order
+
rho => s % rho % array
rho_edge => s % rho_edge % array
rb => grid % rho_base % array
rr => s % rho_p % array
u => s % u % array
+ v => s % v % array
+ kdiff => s % kdiff % array
ru => grid % ru % array
w => s % w % array
rw => grid % rw % array
@@ -1766,6 +1866,7 @@
verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
@@ -1775,6 +1876,10 @@
zz => grid % zz % array
zx => grid % zx % array
+ defc_a => grid % defc_a % array
+ defc_b => grid % defc_b % array
+
+
tend_u => tend % u % array
tend_theta => tend % theta % array
tend_w => tend % w % array
@@ -1846,6 +1951,28 @@
end do
end do
+ delsq_horiz_mixing = .false.
+ if (config_horiz_mixing == "2d_smagorinsky") then
+ do iCell = 1, nCells
+ d_diag(:) = 0.
+ d_off_diag(:) = 0.
+ do iEdge = 1, grid % nEdgesOnCell % array (iCell)
+ do k=1, nVertLevels
+ d_diag(k) = d_diag(k) + defc_a(iEdge,iCell)*u(k,EdgesOnCell(iEdge,iCell)) &
+ - defc_b(iEdge,iCell)*v(k,EdgesOnCell(iEdge,iCell))
+ d_off_diag(k) = d_off_diag(k) + defc_b(iEdge,iCell)*u(k,EdgesOnCell(iEdge,iCell)) &
+ + defc_a(iEdge,iCell)*v(k,EdgesOnCell(iEdge,iCell))
+ end do
+ 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)
+ end do
+ end do
+ delsq_horiz_mixing = .true.
+ else if ( config_horiz_mixing == "2d_fixed") then
+ delsq_horiz_mixing = .true.
+ end if
+
#ifdef LANL_FORMULATION
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
@@ -1942,6 +2069,37 @@
cell2 = cellsOnEdge(2,iEdge)
wduz(1) = 0.
+ if (config_u_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ end do
+
+ else if (config_u_vadv_order == 3) then
+
+ k = 2
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ do k=3,nVertLevels-1
+ wduz(k) = flux3( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)), 1. )
+ end do
+ k = nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+
+ else if (config_u_vadv_order == 4) then
+
+ k = 2
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+ do k=3,nVertLevels-1
+ wduz(k) = flux4( u(k-2,iEdge),u(k-1,iEdge),u(k,iEdge),u(k+1,iEdge),0.5*(rw(k,cell1)+rw(k,cell2)) )
+ end do
+ k = nVertLevels
+ wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
+
+ end if
+
+ wduz(nVertLevels+1) = 0.
+
+ wduz(1) = 0.
do k=2,nVertLevels
wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
end do
@@ -1955,28 +2113,54 @@
!
! horizontal mixing for u
!
- if ( h_mom_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ if (delsq_horiz_mixing) then
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! 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)
- u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
+ if (h_mom_eddy_visc2 > 0.0) then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
+ ! 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)
+ u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
- end do
- end do
- end if
+ else if ( config_horiz_mixing == "2d_smagorinsky") then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
+ ! 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)
+ u_diffusion = rho_edge(k,iEdge)* 0.5*(kdiff(k,cell1)+kdiff(k,cell2)) * u_diffusion
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+ end if
+
+ end if ! delsq_horiz_mixing for u
+
if ( h_mom_eddy_visc4 > 0.0 ) then
allocate(delsq_divergence(nVertLevels, nCells+1))
@@ -2080,7 +2264,7 @@
!
if ( v_mom_eddy_visc2 > 0.0 ) then
- if (mix_full) then
+ if (config_mix_full) then
do iEdge=1,grid % nEdgesSolve
@@ -2144,10 +2328,8 @@
! horizontal advection for w
!
- w_adv_order = 2
+ if (config_w_adv_order == 2) then
- if (w_adv_order == 2) then
-
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2161,7 +2343,7 @@
end if
end do
- else if (w_adv_order == 3) then
+ else if (config_w_adv_order == 3) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2199,7 +2381,7 @@
end if
end do
- else if (w_adv_order == 4) then
+ else if (config_w_adv_order == 4) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2239,9 +2421,11 @@
! Note: we are using quite a bit of the theta code here - could be combined later???
- if ( h_mom_eddy_visc2 > 0.0 ) then
+ if (delsq_horiz_mixing) then
- do iEdge=1,grid % nEdges
+ if (h_mom_eddy_visc2 > 0.0) then
+
+ do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
@@ -2253,8 +2437,27 @@
tend_w(k,cell2) = tend_w(k,cell2) - flux
end do
- end if
- end do
+ end if
+ end do
+ else if (config_horiz_mixing == "2d_smagorinsky") then
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+
+ do k=2,grid % nVertLevels
+! 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)
+ 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
+ end do
+
+ end if
+ end do
+ end if
end if
@@ -2310,14 +2513,40 @@
!
do iCell = 1, nCells
+
wdwz(1) = 0.
- do k=2,nVertLevels
+ if (config_w_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+ end do
+
+ else if (config_w_vadv_order == 3) then
+
+ k = 2
wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
- end do
+ do k=3,nVertLevels-1
+ wdwz(k) = flux3( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)), 1. )
+ end do
+ k = nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+
+ else if (config_w_vadv_order == 4) then
+
+ k = 2
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+ do k=3,nVertLevels-1
+ wdwz(k) = flux4( w(k-2,iCell),w(k-1,iCell),w(k,iCell),w(k+1,iCell),0.5*(rw(k,iCell)+rw(k-1,iCell)) )
+ end do
+ k = nVertLevels
+ wdwz(k) = 0.25*(rw(k,icell)+rw(k-1,iCell))*(w(k,iCell)+w(k-1,iCell))
+
+ end if
+
wdwz(nVertLevels+1) = 0.
+
do k=2,nVertLevels
-
tend_w(k,iCell) = tend_w(k,iCell)/areaCell(iCell) -rdzu(k)*(wdwz(k+1)-wdwz(k)) &
- cqw(k,iCell)*( rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+ gravity* &
@@ -2407,8 +2636,6 @@
! 3rd order stencil
- coef_3rd_order = 0.25
-
if( u(k,iEdge) > 0) then
flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
0.5*(theta(k,cell1) + theta(k,cell2)) &
@@ -2474,22 +2701,42 @@
! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
! but here we can also code in hyperdiffusion if we wish (2nd order at present)
!
- if ( h_theta_eddy_visc2 > 0.0 ) then
+ if (delsq_horiz_mixing) then
+ if ( h_theta_eddy_visc2 > 0.0 ) then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+
+ do k=1,grid % nVertLevels
+ theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(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
+ end do
+
+ end if
+ end do
- do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(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
- end do
+ else if ( ( config_horiz_mixing == "2d_smagorinsky") ) then
- end if
- end do
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ 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(k,cell2) - theta(k,cell1))/dcEdge(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
+ end do
+
+ end if
+ end do
+ end if
end if
@@ -2544,11 +2791,39 @@
! Note: we are also dividing through by the cell area after the horizontal flux divergence
!
do iCell = 1, nCells
+
wdtz(1) = 0.
- do k=2,nVertLevels
+ if (config_theta_vadv_order == 2) then
+
+ do k=2,nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+ end do
+
+ else if (config_theta_vadv_order == 3) then
+
+ k = 2
wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
- end do
+ do k=3,nVertLevels-1
+ wdtz(k) = flux3( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), rw(k,iCell), coef_3rd_order )
+ end do
+ k = nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+
+ else if (config_theta_vadv_order == 4) then
+
+ k = 2
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+ do k=3,nVertLevels-1
+ wdtz(k) = flux4( theta(k-2,iCell),theta(k-1,iCell),theta(k,iCell),theta(k+1,iCell), rw(k,iCell) )
+ end do
+ k = nVertLevels
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta(k,iCell)+fzp(k)*theta(k-1,iCell))
+
+ end if
+
+
wdtz(nVertLevels+1) = 0.
+
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdzw(k)*(wdtz(k+1)-wdtz(k))
!! tend_theta(k,iCell) = tend_theta(k) + rho(k,iCell)*h_diabatic(k,iCell)
@@ -2560,7 +2835,7 @@
!
if ( v_theta_eddy_visc2 > 0.0 ) then
- if (mix_full) then
+ if (config_mix_full) then
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels-1
@@ -3062,6 +3337,6 @@
end do
end do
- end subroutine kessler
+ end subroutine kessler
end module time_integration
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-07-26 16:16:00 UTC (rev 443)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-07-28 21:30:42 UTC (rev 444)
@@ -29,6 +29,7 @@
call init_coupled_diagnostics( block % time_levs(1) % state, mesh)
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh) ! ok for nonhydrostatic model
call initialize_advection_rk(mesh)
+ call initialize_deformation_weights(mesh)
end subroutine mpas_init
</font>
</pre>