<p><b>duda</b> 2012-01-18 10:45:33 -0700 (Wed, 18 Jan 2012)</p><p>BRANCH COMMIT<br>
<br>
Add new horizontal pressure gradient code from Joe.<br>
<br>
<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
M src/core_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/mpas_atm_mpas_core.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry 2012-01-17 22:48:31 UTC (rev 1386)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry 2012-01-18 17:45:33 UTC (rev 1387)
@@ -150,6 +150,8 @@
var persistent real zf3 ( nVertLevelsP1 TWO nEdges ) 0 iro zf3 mesh - -
var persistent real zb ( nVertLevelsP1 TWO nEdges ) 0 iro zb mesh - -
var persistent real zb3 ( nVertLevelsP1 TWO nEdges ) 0 iro zb3 mesh - -
+var persistent real pzm ( nVertLevels nCells ) 0 r pzm mesh - -
+var persistent real pzp ( nVertLevels nCells ) 0 r pzp mesh - -
% coefficients for the vertical tridiagonal solve
% Note: these could be local but...
Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_mpas_core.F 2012-01-17 22:48:31 UTC (rev 1386)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_mpas_core.F 2012-01-18 17:45:33 UTC (rev 1387)
@@ -560,12 +560,14 @@
type (mesh_type), intent(inout) :: mesh
- integer :: iEdge, iCell1, iCell2, k
+ integer :: iEdge, iCell1, iCell2, k, iCell, nz, nz1
real (kind=RKIND) :: d1, d2, d3
- real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl, zgrid
+ real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl, zgrid, pzp, pzm
cpr => mesh % cpr % array
cpl => mesh % cpl % array
+ pzp => mesh % pzp % array
+ pzm => mesh % pzm % array
zgrid => mesh % zgrid % array
!**** coefficient arrays for new pressure gradient calculation
@@ -575,28 +577,70 @@
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))
+! 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))
+ cpr(1,iEdge) = d2/(d2-d1)
+ cpr(2,iEdge) = -d1/(d2-d1)
+ cpr(3,iEdge) = 0.
+
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))
+! 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))
+ cpl(1,iEdge) = d2/(d2-d1)
+ cpl(2,iEdge) = -d1/(d2-d1)
+ cpl(3,iEdge) = 0.
+
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)
+ else
+
+! Coefficients for computing vertical pressure gradient dp/dz
+! dp/dz (k,iCell) = pzp(k,iCell) * (p(k+1,iCell) - p(k,iCell)) +pzm(k,iCell) * (p(k,iCell) - p(k-1,iCell))
+
+ nz1 = mesh % nVertLevels
+ nz = nz1 + 1
+
+ do iCell=1, mesh % nCells
+
+ d1 = zgrid(3,iCell)-zgrid(1,iCell)
+ d2 = zgrid(4,iCell)-zgrid(2,iCell)
+ d3 = d1+d2
+ pzm(1,iCell) = 2.*d3/(d1*d2)
+ pzp(1,iCell) = -2.*d1/(d2*d3)
+
+ do k=2,nz1-1
+ pzp(k,iCell) = 2.*(zgrid(k+1,iCell)-zgrid(k-1,iCell))/ &
+ & ((zgrid(k+2,iCell)-zgrid(k ,iCell))* &
+ & (zgrid(k+2,iCell)-zgrid(k ,iCell) &
+ & +zgrid(k+1,iCell)-zgrid(k-1,iCell)))
+ pzm(k,iCell) = 2.*(zgrid(k+2,iCell)-zgrid(k ,iCell))/ &
+ & ((zgrid(k+1,iCell)-zgrid(k-1,iCell))* &
+ & (zgrid(k+2,iCell)-zgrid(k ,iCell) &
+ & +zgrid(k+1,iCell)-zgrid(k-1,iCell)))
+ end do
+
+ pzp(nz1,iCell) = 0.
+ pzm(nz1,iCell) = 2./(zgrid(nz,iCell)-zgrid(nz1-1,iCell))
+
+ end do
+
end if
end subroutine atm_compute_pgf_coefs
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-01-17 22:48:31 UTC (rev 1386)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-01-18 17:45:33 UTC (rev 1387)
@@ -749,7 +749,7 @@
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), dimension(:,:), pointer :: cpr, cpl, pzp, pzm
real (kind=RKIND) :: smdiv, c2, rcv
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
@@ -797,6 +797,9 @@
gamma_tri => diag % gamma_tri % array
dss => grid % dss % array
+ pzp => grid % pzp % array
+ pzm => grid % pzm % array
+
tend_ru => tend % u % array
tend_rho => tend % rho_zz % array
tend_rt => tend % theta_m % array
@@ -883,23 +886,46 @@
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)))
+! 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)))
+
+ dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
+ *(pzm(k,cell2)*(zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
+ -zz(k ,cell2)*rtheta_pp_old(k ,cell2)) &
+ +pzm(k,cell1)*(zz(k+1,cell1)*rtheta_pp_old(k+1,cell1) &
+ -zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
+ +pzp(k,cell2)*(zz(k+2,cell2)*rtheta_pp_old(k+2,cell2) &
+ -zz(k ,cell2)*rtheta_pp_old(k ,cell2)) &
+ +pzp(k,cell1)*(zz(k+2,cell1)*rtheta_pp_old(k+2,cell1) &
+ -zz(k ,cell1)*rtheta_pp_old(k ,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)))
+! 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)))
+ dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
+ *(pzp(k,cell2)*(zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
+ -zz(k ,cell2)*rtheta_pp_old(k ,cell2)) &
+ +pzm(k,cell2)*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
+ -zz(k-1,cell2)*rtheta_pp_old(k-1,cell2)) &
+ +pzp(k,cell1)*(zz(k+1,cell1)*rtheta_pp_old(k+1,cell1) &
+ -zz(k ,cell1)*rtheta_pp_old(k ,cell1)) &
+ +pzm(k,cell1)*(zz(k ,cell1)*rtheta_pp_old(k ,cell1) &
+ -zz(k-1,cell1)*rtheta_pp_old(k-1,cell1)))
end do
- dpzx(nVertLevels + 1) = 0.
+! 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 = (rtheta_pp_old(k,cell2)-rtheta_pp_old(k,cell1))/dcEdge(iEdge) &
+! - rdzw(k)*(dpzx(k+1)-dpzx(k))
+ pgrad = ((rtheta_pp_old(k,cell2)*zz(k,cell2) &
+ -rtheta_pp_old(k,cell1)*zz(k,cell1))/dcEdge(iEdge) &
+ -dpzx(k))/(.5*(zz(k,cell2)+zz(k,cell1)))
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))
@@ -1841,7 +1867,7 @@
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
+ real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl, pzp, pzm
integer :: kr, kl
real (kind=RKIND), allocatable, dimension(:,:) :: rv, divergence_ru, qtot
@@ -1907,7 +1933,10 @@
pressure_b => diag % pressure_base % array
h_divergence => diag % h_divergence % array
+ pzp => grid % pzp % array
+ pzm => grid % pzm % array
+
weightsOnEdge => grid % weightsOnEdge % array
cellsOnEdge => grid % cellsOnEdge % array
verticesOnEdge => grid % verticesOnEdge % array
@@ -2046,16 +2075,64 @@
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.
+ if(newpx) then
+ 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
+
+ else
+ 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)))
+
+ dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
+ *(pzm(k,cell2)*(pp(k+1,cell2)-pp(k,cell2)) &
+ +pzm(k,cell1)*(pp(k+1,cell1)-pp(k,cell1)) &
+ +pzp(k,cell2)*(pp(k+2,cell2)-pp(k,cell2)) &
+ +pzp(k,cell1)*(pp(k+2,cell1)-pp(k,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)))
+
+ dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
+ *(pzp(k,cell2)*(pp(k+1,cell2)-pp(k ,cell2)) &
+ +pzm(k,cell2)*(pp(k ,cell2)-pp(k-1,cell2)) &
+ +pzp(k,cell1)*(pp(k+1,cell1)-pp(k ,cell1)) &
+ +pzm(k,cell1)*(pp(k ,cell1)-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)) )
+
+ tend_u(k,iEdge) = - cqu(k,iEdge)*((pp(k,cell2)-pp(k,cell1))/dcEdge(iEdge) &
+ - dpzx(k) ) / (.5*(zz(k,cell2)+zz(k,cell1)))
+ end do
+
+ end if
+
wduz(1) = 0.
k = 2
wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
@@ -2068,8 +2145,8 @@
wduz(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)) )
+! 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)) )
tend_u(k,iEdge) = tend_u(k,iEdge) - rdzw(k)*(wduz(k+1)-wduz(k))
end do
</font>
</pre>