<p><b>duda</b> 2011-08-26 15:10:30 -0600 (Fri, 26 Aug 2011)</p><p>Add code for new horizontal pressure gradient calculation. Whether the new code <br>
or the old code is used can be selected at run-time via the namelist switch <br>
config_newpx (.true. -> new code, .false. -> old code).<br>
<br>
Also, add the &damping namelist record for absorbing layer from r951 to all <br>
non-hydrostatic namelists, since many compilers create code that stops when <br>
a namelist record isn't found.<br>
<br>
<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_squall<br>
M namelist.input.nhyd_atmos<br>
M namelist.input.nhyd_atmos_jw<br>
M namelist.input.nhyd_atmos_mtn_wave<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/namelist.input.nhyd_atmos
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/namelist.input.nhyd_atmos 2011-08-26 21:10:30 UTC (rev 957)
@@ -19,8 +19,14 @@
config_epssm = 0.1
config_smdiv = 0.1
config_h_ScaleWithMesh = .false.
+ config_newpx = .true.
/
+&damping
+ config_zd = 22000.0
+ config_xnutr = 0.0
+/
+
&dimensions
config_nvertlevels = 26
/
Modified: branches/atmos_physics/namelist.input.nhyd_atmos_jw
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos_jw 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/namelist.input.nhyd_atmos_jw 2011-08-26 21:10:30 UTC (rev 957)
@@ -26,8 +26,14 @@
config_epssm = 0.1
config_smdiv = 0.1
config_h_ScaleWithMesh = .false.
+ config_newpx = .true.
/
+&damping
+ config_zd = 22000.0
+ config_xnutr = 0.0
+/
+
&dimensions
config_nvertlevels = 26
/
Modified: branches/atmos_physics/namelist.input.nhyd_atmos_mtn_wave
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos_mtn_wave 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/namelist.input.nhyd_atmos_mtn_wave 2011-08-26 21:10:30 UTC (rev 957)
@@ -20,8 +20,15 @@
config_horiz_mixing = '2d_fixed'
config_coef_3rd_order = 0.25
config_epssm = 0.2
+ config_newpx = .true.
/
+
+&damping
+ config_zd = 22000.0
+ config_xnutr = 0.0
+/
+
&dimensions
config_nvertlevels = 26
/
Modified: branches/atmos_physics/namelist.input.nhyd_atmos_squall
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_atmos_squall 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/namelist.input.nhyd_atmos_squall 2011-08-26 21:10:30 UTC (rev 957)
@@ -15,8 +15,14 @@
config_scalar_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
+ config_newpx = .true.
/
+&damping
+ config_zd = 22000.0
+ config_xnutr = 0.0
+/
+
&dimensions
config_nvertlevels = 26
/
Modified: branches/atmos_physics/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-08-26 21:10:30 UTC (rev 957)
@@ -31,6 +31,7 @@
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 logical nhyd_model config_newpx true
namelist real nhyd_model config_apvm_upwinding 0.5
namelist logical nhyd_model config_h_ScaleWithMesh false
namelist real damping config_zd 22000.0
@@ -127,6 +128,9 @@
var persistent real cf2 ( ) 0 iro cf2 mesh - -
var persistent real cf3 ( ) 0 iro cf3 mesh - -
+var persistent real cpr ( THREE nEdges ) 0 ro cpr mesh - -
+var persistent real cpl ( THREE nEdges ) 0 ro cpl mesh - -
+
# description of the vertical grid structure
var persistent real hx ( nVertLevelsP1 nCells ) 0 iro hx mesh - -
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-08-26 21:10:30 UTC (rev 957)
@@ -124,6 +124,8 @@
call compute_damping_coefs(mesh)
+ call compute_pgf_coefs(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)), &
@@ -374,4 +376,55 @@
end subroutine compute_damping_coefs
+
+ subroutine compute_pgf_coefs(mesh)
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer :: iEdge, iCell1, iCell2, k
+ real (kind=RKIND) :: d1, d2, d3
+ real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl, zgrid
+
+ cpr => mesh % cpr % array
+ cpl => mesh % cpl % array
+ zgrid => mesh % zgrid % array
+
+!**** coefficient arrays for new pressure gradient calculation
+
+ cpr(:,:) = 0.0
+ cpl(:,:) = 0.0
+
+ if (config_newpx) then
+ do iEdge=1,mesh%nEdges
+ iCell1 = mesh % cellsOnEdge % array(1,iEdge)
+ iCell2 = mesh % cellsOnEdge % array(2,iEdge)
+
+ d1 = .25*(zgrid(1,iCell2)+zgrid(2,iCell2)-zgrid(1,iCell1)-zgrid(2,iCell1))
+ d2 = d1+.5*(zgrid(3,iCell2)-zgrid(1,iCell2))
+ d3 = d2+.5*(zgrid(4,iCell2)-zgrid(2,iCell2))
+ cpr(1,iEdge) = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ cpr(2,iEdge) = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ cpr(3,iEdge) = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ d1 = .25*(zgrid(1,iCell1)+zgrid(2,iCell1)-zgrid(1,iCell2)-zgrid(2,iCell2))
+ d2 = d1+.5*(zgrid(3,iCell1)-zgrid(1,iCell1))
+ d3 = d2+.5*(zgrid(4,iCell1)-zgrid(2,iCell1))
+ cpl(1,iEdge) = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ cpl(2,iEdge) = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ cpl(3,iEdge) = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ end do
+! write(6,*) 'cpr1 = ',cpr(1,1),' cpl1 = ',cpl(1,1)
+! write(6,*) 'cpr2 = ',cpr(2,1),' cpl2 = ',cpl(2,1)
+! write(6,*) 'cpr3 = ',cpr(3,1),' cpl3 = ',cpl(3,1)
+
+ end if
+
+ end subroutine compute_pgf_coefs
+
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-08-26 17:47:32 UTC (rev 956)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-08-26 21:10:30 UTC (rev 957)
@@ -735,6 +735,8 @@
zgrid, cofwr, cofwz, w, h_divergence
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, rdzw, dcEdge, AreaCell, cofrz, dvEdge
+ real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl
+
real (kind=RKIND) :: smdiv, c2, rcv
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: dpzx
@@ -743,7 +745,8 @@
integer :: cell1, cell2, iEdge, iCell, k
real (kind=RKIND) :: pgrad, flux1, flux2, flux, resm, epssm
- real (kind=RKIND) :: cf1, cf2, cf3
+ real (kind=RKIND) :: cf1, cf2, cf3, pr, pl
+ integer :: kr, kl
integer :: nEdges, nCells, nCellsSolve, nVertLevels
@@ -752,6 +755,7 @@
logical, parameter :: debug1 = .false.
real (kind=RKIND) :: wmax
integer :: iwmax, kwmax
+ logical :: newpx
!--
@@ -805,6 +809,10 @@
cf2 = grid % cf2 % scalar
cf3 = grid % cf3 % scalar
+ cpr => grid % cpr % array
+ cpl => grid % cpl % array
+ newpx = config_newpx
+
epssm = config_epssm
smdiv = config_smdiv
@@ -822,54 +830,90 @@
do iEdge = 1, nEdges
- cell1 = grid % cellsOnEdge % array (1,iEdge)
- cell2 = grid % cellsOnEdge % array (2,iEdge)
- ! update edge for block-owned cells
- if (cell1 <= grid % nCellsSolve .or. cell2 <= grid % nCellsSolve ) then
+ cell1 = grid % cellsOnEdge % array (1,iEdge)
+ cell2 = grid % cellsOnEdge % array (2,iEdge)
+ ! update edge for block-owned cells
+ if (cell1 <= grid % nCellsSolve .or. cell2 <= grid % nCellsSolve ) then
- k = 1
- dpzx(k) = .5*zx(k,iEdge)*(cf1*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
- +zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
- +cf2*(zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
- +zz(k+1,cell1)*rtheta_pp_old(k+1,cell1)) &
- +cf3*(zz(k+2,cell2)*rtheta_pp_old(k+2,cell2) &
- +zz(k+2,cell1)*rtheta_pp_old(k+2,cell1)))
- do k=2,grid % nVertLevels
- dpzx(k)=.5*zx(k,iEdge)*(fzm(k)*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
- +zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
- +fzp(k)*(zz(k-1,cell2)*rtheta_pp_old(k-1,cell2) &
- +zz(k-1,cell1)*rtheta_pp_old(k-1,cell1)))
- end do
- dpzx(nVertLevels + 1) = 0.
+ if (newpx) then
- do k=1,nVertLevels
- pgrad = (rtheta_pp_old(k,cell2)-rtheta_pp_old(k,cell1))/dcEdge(iEdge) &
- - 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.05/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
+ k = 1
+ pr = cpr(k ,iEdge)*zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
+ + cpr(k+1,iEdge)*zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
+ + cpr(k+2,iEdge)*zz(k+2,cell2)*rtheta_pp_old(k+2,cell2)
- ru_p(k,iEdge) = ru_p(k,iEdge) + du(k)
+ pl = cpl(k ,iEdge)*zz(k ,cell1)*rtheta_pp_old(k ,cell1) &
+ + cpl(k+1,iEdge)*zz(k+1,cell1)*rtheta_pp_old(k+1,cell1) &
+ + cpl(k+2,iEdge)*zz(k+2,cell1)*rtheta_pp_old(k+2,cell1)
+ pgrad = 2./(zz(k,cell1)+zz(k,cell2))*(pr-pl)/dcEdge(iEdge)
+ pgrad = 0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad
+ du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
- if(debug) then
- if(iEdge == 3750) then
- write(0,*) ' k, pgrad, tend_ru ',k,pgrad,tend_ru(k,3750)
- end if
+ do k=2,nVertLevels
+
+ kr = min(nVertLevels,k+ nint(.5-sign(.5,zx(k,iEdge)+zx(k+1,iEdge))))
+ kl = min(nVertLevels,2*k+1-kr)
+ pr = zz(k,cell2)*rtheta_pp_old(k ,cell2)+.5*(zgrid(k ,cell1) +zgrid(k +1,cell1) &
+ -zgrid(k ,cell2) -zgrid(k +1,cell2)) &
+ /(zgrid(kr+1,cell2) -zgrid(kr-1,cell2)) &
+ *(zz(kr,cell2)*rtheta_pp_old(kr,cell2)-zz(kr-1,cell2)*rtheta_pp_old(kr-1,cell2))
+ pl = zz(k,cell1)*rtheta_pp_old(k ,cell1)+.5*(zgrid(k ,cell2) +zgrid(k +1,cell2) &
+ -zgrid(k ,cell1) -zgrid(k +1,cell1)) &
+ /(zgrid(kl+1,cell1) -zgrid(kl-1,cell1)) &
+ *(zz(kl,cell1)*rtheta_pp_old(kl,cell1)-zz(kl-1,cell1)*rtheta_pp_old(kl-1,cell1))
+ pgrad = 2./(zz(k,cell1)+zz(k,cell2))*(pr-pl)/dcEdge(iEdge)
+ pgrad = 0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad
+ du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
+ end do
+
+ else
+
+ k = 1
+ dpzx(k) = .5*zx(k,iEdge)*(cf1*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
+ +zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
+ +cf2*(zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
+ +zz(k+1,cell1)*rtheta_pp_old(k+1,cell1)) &
+ +cf3*(zz(k+2,cell2)*rtheta_pp_old(k+2,cell2) &
+ +zz(k+2,cell1)*rtheta_pp_old(k+2,cell1)))
+ do k=2,grid % nVertLevels
+ dpzx(k)=.5*zx(k,iEdge)*(fzm(k)*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
+ +zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
+ +fzp(k)*(zz(k-1,cell2)*rtheta_pp_old(k-1,cell2) &
+ +zz(k-1,cell1)*rtheta_pp_old(k-1,cell1)))
+ end do
+ dpzx(nVertLevels + 1) = 0.
+
+ do k=1,nVertLevels
+ pgrad = (rtheta_pp_old(k,cell2)-rtheta_pp_old(k,cell1))/dcEdge(iEdge) &
+ - 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.05/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
+ end do
end if
-! need to add horizontal fluxes into density update, rtheta update and w update
+ do k=1,nVertLevels
+ ru_p(k,iEdge) = ru_p(k,iEdge) + du(k)
- flux = dts*dvEdge(iEdge)*ru_p(k,iEdge)
- rs(k,cell1) = rs(k,cell1)-flux/AreaCell(cell1)
- rs(k,cell2) = rs(k,cell2)+flux/AreaCell(cell2)
+ if(debug) then
+ if(iEdge == 3750) then
+ write(0,*) ' k, pgrad, tend_ru ',k,pgrad,tend_ru(k,3750)
+ end if
+ end if
- flux = flux*0.5*(theta_m(k,cell2)+theta_m(k,cell1))
- ts(k,cell1) = ts(k,cell1)-flux/AreaCell(cell1)
- ts(k,cell2) = ts(k,cell2)+flux/AreaCell(cell2)
+! need to add horizontal fluxes into density update, rtheta update and w update
- ruAvg(k,iEdge) = ruAvg(k,iEdge) + ru_p(k,iEdge)
+ flux = dts*dvEdge(iEdge)*ru_p(k,iEdge)
+ rs(k,cell1) = rs(k,cell1)-flux/AreaCell(cell1)
+ rs(k,cell2) = rs(k,cell2)+flux/AreaCell(cell2)
+
+ flux = flux*0.5*(theta_m(k,cell2)+theta_m(k,cell1))
+ ts(k,cell1) = ts(k,cell1)-flux/AreaCell(cell1)
+ ts(k,cell2) = ts(k,cell2)+flux/AreaCell(cell2)
+
+ ruAvg(k,iEdge) = ruAvg(k,iEdge) + ru_p(k,iEdge)
- end do
+ end do
end if ! end test for block-owned cells
@@ -1956,11 +2000,14 @@
real (kind=RKIND), dimension(:), pointer :: rdzu, rdzw, fzm, fzp, qv_init
real (kind=RKIND), dimension(:,:), pointer :: t_init
+ real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl
+ integer :: kr, kl
+
real (kind=RKIND), allocatable, dimension(:,:) :: rv, divergence_ru, qtot
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_theta, delsq_divergence
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
- real (kind=RKIND) :: cf1, cf2, cf3
+ real (kind=RKIND) :: cf1, cf2, cf3, pr, pl
! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
@@ -1974,7 +2021,7 @@
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
+ logical :: delsq_horiz_mixing, newpx
real (kind=RKIND) :: coef_3rd_order
@@ -2062,6 +2109,10 @@
cqw => diag % cqw % array
cqu => diag % cqu % array
+ cpr => grid % cpr % array
+ cpl => grid % cpl % array
+ newpx = config_newpx
+
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -2141,23 +2192,64 @@
delsq_horiz_mixing = .true.
end if
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+!**** horizontal pressure gradient ***************************************************
- ! horizontal pressure gradient, nonlinear Coriolis term and ke gradient
+ if (newpx) then
- k = 1
- dpzx(k) = .5*zx(k,iEdge)*(cf1*(pp(k ,cell2)+pp(k ,cell1)) &
- +cf2*(pp(k+1,cell2)+pp(k+1,cell1)) &
- +cf3*(pp(k+2,cell2)+pp(k+2,cell1)))
- do k = 2, nVertLevels
- dpzx(k) = .5*zx(k,iEdge)*(fzm(k)*(pp(k ,cell2)+pp(k ,cell1)) &
- +fzp(k)*(pp(k-1,cell2)+pp(k-1,cell1)))
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ k = 1
+ pr = cpr(k,iEdge)*pp(k,cell2)+cpr(k+1,iEdge)*pp(k+1,cell2)+cpr(k+2,iEdge)*pp(k+2,cell2)
+ pl = cpl(k,iEdge)*pp(k,cell1)+cpl(k+1,iEdge)*pp(k+1,cell1)+cpl(k+2,iEdge)*pp(k+2,cell1)
+ tend_u(k,iEdge) = - cqu(k,iEdge)*2./(zz(k,cell1)+zz(k,cell2))*(pr-pl)/dcEdge(iEdge)
+
+ do k=2,nVertLevels
+
+ kr = min(nVertLevels,k+ nint(.5-sign(.5,zx(k,iEdge)+zx(k+1,iEdge))))
+ kl = min(nVertLevels,2*k+1-kr)
+
+ pr = pp(k,cell2)+.5*(zgrid(k ,cell1)+zgrid(k +1,cell1)-zgrid(k ,cell2)-zgrid(k +1,cell2)) &
+ /(zgrid(kr+1,cell2)-zgrid(kr-1,cell2))*( pp(kr,cell2)-pp (kr-1,cell2))
+ pl = pp(k,cell1)+.5*(zgrid(k ,cell2)+zgrid(k +1,cell2)-zgrid(k ,cell1)-zgrid(k +1,cell1)) &
+ /(zgrid(kl+1,cell1)-zgrid(kl-1,cell1))*( pp(kl,cell1)-pp (kl-1,cell1))
+ tend_u(k,iEdge) = - cqu(k,iEdge)*2./(zz(k,cell1)+zz(k,cell2))*(pr-pl)/dcEdge(iEdge)
+
+ end do
end do
- dpzx(nVertLevels+1) = 0.
+ else
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ k = 1
+ dpzx(k) = .5*zx(k,iEdge)*(cf1*(pp(k ,cell2)+pp(k ,cell1)) &
+ +cf2*(pp(k+1,cell2)+pp(k+1,cell1)) &
+ +cf3*(pp(k+2,cell2)+pp(k+2,cell1)))
+ do k=2,nVertLevels
+ dpzx(k) = .5*zx(k,iEdge)*(fzm(k)*(pp(k ,cell2)+pp(k ,cell1)) &
+ +fzp(k)*(pp(k-1,cell2)+pp(k-1,cell1)))
+ end do
+ dpzx(nVertLevels+1) = 0.
+
+ do k=1,nVertLevels
+ tend_u(k,iEdge) = - cqu(k,iEdge)*( (pp(k,cell2)/zz(k,cell2) - pp(k,cell1)/zz(k,cell1)) &
+ / dcEdge(iEdge) - rdzw(k)*(dpzx(k+1)-dpzx(k)) )
+ end do
+ end do
+
+ end if
+
+!**** nonlinear Coriolis term and ke gradient ****************************************
+
+
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
do k=1,nVertLevels
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
@@ -2165,10 +2257,9 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * rho_edge(k,eoe)
end do
- tend_u(k,iEdge) = rho_edge(k,iEdge)* (q - (ke(k,cell2) - ke(k,cell1)) / dcEdge(iEdge)) &
- - u(k,iEdge)*0.5*(divergence_ru(k,cell1)+divergence_ru(k,cell2)) &
- - cqu(k,iEdge)*( (pp(k,cell2)/zz(k,cell2) - pp(k,cell1)/zz(k,cell1)) / dcEdge(iEdge) &
- -rdzw(k)*(dpzx(k+1)-dpzx(k)) )
+ tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge)* (q - (ke(k,cell2) - ke(k,cell1)) &
+ / dcEdge(iEdge)) &
+ - u(k,iEdge)*0.5*(divergence_ru(k,cell1)+divergence_ru(k,cell2))
end do
end do
</font>
</pre>