<p><b>duda</b> 2011-08-24 17:29:04 -0600 (Wed, 24 Aug 2011)</p><p>BRANCH COMMIT<br>
<br>
Modify the non-hydrostatic core and initialization so that the fields<br>
named 'theta' and 'rho' in the input and output files are potential<br>
temperature and dry density; theta and rho have been added to the<br>
Registry file, and are read from input and written to output. The<br>
prognostic fields which represent modified moist potential temperature<br>
and dry density divided by zeta_z have been renamed 'theta_m' and<br>
'rho_zz' throughout the code, and appear in the restart files as theta_m<br>
and rho_zz.<br>
<br>
NB: Initial conditions must provide 'theta' and 'rho', rather than<br>
'theta_m' and 'rho_zz', since the latter two are re-computed on model<br>
cold starts.<br>
<br>
The correspondence between fields in input and output files before and<br>
after this change is summarized below.<br>
<br>
before after field contents<br>
--------------+---------------+-----------------<br>
theta theta_m modified moist potential temperature<br>
rho rho_zz dry density / zeta_z<br>
[not avail.] theta potential temperature<br>
[not avail.] rho dry density<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M src/core_physics/module_physics_todynamics.F<br>
M src/core_physics/module_physics_interface_nhyd.F<br>
M src/core_nhyd_atmos/module_mpas_core.F<br>
M src/core_nhyd_atmos/module_test_cases.F<br>
M src/core_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_time_integration.F<br>
M graphics/ncl/cells_nhyd_sphere.ncl<br>
M graphics/ncl/cells_nhyd_sph1.ncl<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/graphics/ncl/cells_nhyd_sph1.ncl
===================================================================
--- branches/atmos_physics/graphics/ncl/cells_nhyd_sph1.ncl 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/graphics/ncl/cells_nhyd_sph1.ncl 2011-08-24 23:29:04 UTC (rev 953)
@@ -61,6 +61,7 @@
latEdge = f->latEdge(:) * r2d
verticesOnCell = f->verticesOnCell(:,:)
alpha = f->angleEdge(:)
+ zz = f->zz(:,:)
res = True
res@gsnMaximize = True
@@ -142,8 +143,8 @@
pthird = f->pressure_p(t,:,2)+f->pressure_base(t,:,2)
; fld = (cf1*pfirst + cf2*psecond + cf3*pthird)/100.
- rhofirst = f->rho(t,:,0)
- rhosecond = f->rho(t,:,1)
+ rhofirst = f->rho(t,:,0) / zz(:,0)
+ rhosecond = f->rho(t,:,1) / zz(:,1)
qvfirst = f->qv(t,:,0)
qvsecond = f->qv(t,:,1)
rdzw = f->rdzw
Modified: branches/atmos_physics/graphics/ncl/cells_nhyd_sphere.ncl
===================================================================
--- branches/atmos_physics/graphics/ncl/cells_nhyd_sphere.ncl 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/graphics/ncl/cells_nhyd_sphere.ncl 2011-08-24 23:29:04 UTC (rev 953)
@@ -61,6 +61,7 @@
latEdge = f->latEdge(:) * r2d
verticesOnCell = f->verticesOnCell(:,:)
alpha = f->angleEdge(:)
+ zz = f->zz(:,:)
res = True
res@gsnMaximize = True
@@ -149,8 +150,8 @@
pthird = f->pressure_p(t,:,2)+f->pressure_base(t,:,2)
fld = (cf1*pfirst + cf2*psecond + cf3*pthird)/100.
- rhofirst = f->rho(t,:,0)
- rhosecond = f->rho(t,:,1)
+ rhofirst = f->rho(t,:,0) / zz(:,0)
+ rhosecond = f->rho(t,:,1) / zz(:,1)
qvfirst = f->qv(t,:,0)
qvsecond = f->qv(t,:,1)
rdzw = f->rdzw
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-08-24 23:29:04 UTC (rev 953)
@@ -170,8 +170,8 @@
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 o u state - -
var persistent real w ( nVertLevelsP1 nCells Time ) 2 o w state - -
-var persistent real rho ( nVertLevels nCells Time ) 2 o rho state - -
-var persistent real theta ( nVertLevels nCells Time ) 2 o theta state - -
+var persistent real rho_zz ( nVertLevels nCells Time ) 2 o rho_zz state - -
+var persistent real theta_m ( nVertLevels nCells Time ) 2 o theta_m state - -
var persistent real qv ( nVertLevels nCells Time ) 2 o qv state scalars moist
var persistent real qc ( nVertLevels nCells Time ) 2 o qc state scalars moist
var persistent real qr ( nVertLevels nCells Time ) 2 o qr state scalars moist
@@ -184,6 +184,8 @@
var persistent real qv_init ( nVertLevels ) 0 io qv_init mesh - -
# Diagnostic fields: only written to output
+var persistent real rho ( nVertLevels nCells Time ) 1 o rho diag - -
+var persistent real theta ( nVertLevels nCells Time ) 1 o theta diag - -
var persistent real v ( nVertLevels nEdges Time ) 1 o v diag - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 1 o uReconstructX diag - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 1 o uReconstructY diag - -
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -149,7 +149,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
- real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
+ real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt
real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3, zb, zb3
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
@@ -248,9 +248,9 @@
ppb => diag % pressure_base % array
pp => diag % pressure_p % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
scalars => state % scalars % array
@@ -594,7 +594,7 @@
p (k,i) = ((ppb(k,i)+pp(k,i))/p0)**(rgas/cp)
t (k,i) = tt(k)/p(k,i)
rt (k,i) = t(k,i)*rr(k,i)+rb(k,i)*(t(k,i)-tb(k,i))
- rho (k,i) = rb(k,i) + rr(k,i)
+ rho_zz (k,i) = rb(k,i) + rr(k,i)
end do
if(i == 1) then
@@ -655,7 +655,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- diag % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+ diag % ru % array (k,iEdge) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*state % u % array (k,iEdge)
end do
end if
@@ -761,11 +761,11 @@
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
do k=2,grid%nVertLevels
state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
end do
end do
@@ -795,6 +795,14 @@
write(0,*) ' i, psurf, lat ',i,psurf,grid%latCell%array(i)*180./3.1415828
enddo
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(state % index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_jw
subroutine calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
@@ -864,7 +872,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
!This is temporary variable here. It just need when calculate tangential velocity v.
@@ -943,11 +951,11 @@
p => diag % exner % array
cqw => diag % cqw % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
pp => diag % pressure_p % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
u => state % u % array
ru => diag % ru % array
@@ -1313,7 +1321,7 @@
do i=1,grid % ncells
do k=1,nz1
grid % t_init % array(k,i) = t_init_1d(k)
- rho(k,i) = rb(k,i)+rr(k,i)
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
end do
end do
@@ -1322,7 +1330,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- ru (k,i) = 0.5*(rho(k,cell1)+rho(k,cell2))*u(k,i)
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
end do
end if
end do
@@ -1371,6 +1379,14 @@
! write(0,'(i2,3(2x,f14.10)') k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
! end do
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(state % index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_squall_line
@@ -1393,7 +1409,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zf, zf3, zb, zb3
!This is temporary variable here. It just need when calculate tangential velocity v.
@@ -1486,11 +1502,11 @@
p => diag % exner % array
cqw => diag % cqw % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
pp => diag % pressure_p % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
u => state % u % array
ru => diag % ru % array
@@ -1833,7 +1849,7 @@
do i=1,grid % ncells
do k=1,nz1
- rho(k,i) = rb(k,i)+rr(k,i)
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
end do
do k=1,nz1
@@ -1846,7 +1862,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- ru (k,i) = 0.5*(rho(k,cell1)+rho(k,cell2))*u(k,i)
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
end do
end if
end do
@@ -1937,11 +1953,11 @@
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
do k=2,grid%nVertLevels
state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
end do
end do
@@ -1973,6 +1989,14 @@
! write(0,*) ' k,u_init, t_init, qv_init ',k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
! end do
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(state % index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_mtn_wave
@@ -2012,7 +2036,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:), pointer :: vert_level, latPoints, lonPoints, ter
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
- real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
+ real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt
real (kind=RKIND), dimension(:), pointer :: destField1d
real (kind=RKIND), dimension(:,:), pointer :: destField2d
real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3, zb, zb3
@@ -2133,9 +2157,9 @@
ppb => diag % pressure_base % array
pp => diag % pressure_p % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
scalars => state % scalars % array
@@ -3695,7 +3719,7 @@
call quicksort(config_nfglevels, sorted_arr)
do k=1,grid%nVertLevels
target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
- state % theta % array(k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1)
+ state % theta_m % array(k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1)
end do
@@ -3816,12 +3840,13 @@
!
! Diagnose fields needed in initial conditions file (u, w, rho, theta, scalars)
+ ! NB: At this point, "rho_zz" is simple dry density, and "theta_m" is regular potential temperature
!
do iCell=1,grid%nCells
do k=1,grid%nVertLevels
! QV
- es = 6.112 * exp((17.27*(state % theta % array(k,iCell) - 273.16))/(state % theta % array(k,iCell) - 35.86))
+ es = 6.112 * exp((17.27*(state % theta_m % array(k,iCell) - 273.16))/(state % theta_m % array(k,iCell) - 35.86))
rs = 0.622 * es / (diag % pressure % array(k,iCell) - es)
scalars(state % index_qv,k,iCell) = rs * scalars(state % index_qv,k,iCell)
@@ -3832,9 +3857,9 @@
! t(k,iCell) = t(k,iCell) / p(k,iCell)
t(k,iCell) = t(k,iCell) * (p0 / diag % pressure % array(k,iCell)) ** (rgas / cp)
- ! RHO
- rho(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell) * t(k,iCell))
- rho(k,iCell) = rho(k,iCell) / (1.0 + scalars(state % index_qv,k,iCell))
+ ! RHO_ZZ
+ rho_zz(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell) * t(k,iCell))
+ rho_zz(k,iCell) = rho_zz(k,iCell) / (1.0 + scalars(state % index_qv,k,iCell))
end do
end do
@@ -3860,14 +3885,14 @@
do iCell=1,grid%nCells
do k=1,grid%nVertLevels
pp(k,iCell) = diag % pressure % array(k,iCell) - ppb(k,iCell)
- rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
+ rr(k,iCell) = rho_zz(k,iCell) - rb(k,iCell)
end do
end do
do iCell=1,grid%nCells
k = 1
- rho(k,iCell) = ((diag % pressure % array(k,iCell) / p0)**(cv / cp)) * (p0 / rgas) / (t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
- rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
+ rho_zz(k,iCell) = ((diag % pressure % array(k,iCell) / p0)**(cv / cp)) * (p0 / rgas) / (t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
+ rr(k,iCell) = rho_zz(k,iCell) - rb(k,iCell)
do k=2,grid % nVertLevels
it = 0
@@ -3877,11 +3902,11 @@
p_check = pp(k,iCell)
dz = (zgrid(k,iCell) - zgrid(k-1,iCell))
pp(k,iCell) = pp(k-1,iCell) - 0.5 * (rr(k,iCell) + rr(k-1,iCell))*gravity*dz &
- - 0.5 * (rho(k,iCell)*scalars(state % index_qv,k,iCell) + rho(k-1,iCell)*scalars(state % index_qv,k-1,iCell))*gravity*dz
+ - 0.5 * (rho_zz(k,iCell)*scalars(state % index_qv,k,iCell) + rho_zz(k-1,iCell)*scalars(state % index_qv,k-1,iCell))*gravity*dz
diag % pressure % array(k,iCell) = pp(k,iCell) + ppb(k,iCell)
p(k,iCell) = (diag % pressure % array(k,iCell) / p0) ** (rgas / cp)
- rho(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell)*t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
- rr(k,iCell) = rho(k,iCell) - rb(k,iCell)
+ rho_zz(k,iCell) = diag % pressure % array(k,iCell) / rgas / (p(k,iCell)*t(k,iCell)*(1.0 + 1.61*scalars(state % index_qv,k,iCell)))
+ rr(k,iCell) = rho_zz(k,iCell) - rb(k,iCell)
p_check = abs(p_check - pp(k,iCell))
@@ -3894,14 +3919,14 @@
do iCell=1,grid%nCells
do k=1,grid%nVertLevels
t(k,iCell) = t(k,iCell) * (1.0 + 1.61*scalars(state % index_qv,k,iCell))
- rho(k,iCell) = rho(k,iCell) / zz(k,iCell)
+ rho_zz(k,iCell) = rho_zz(k,iCell) / zz(k,iCell)
rb(k,iCell) = rb(k,iCell) / zz(k,iCell)
end do
end do
do iEdge=1,grid%nEdges
do k=1,grid%nVertLevels
- diag % ru % array(k,iEdge) = state % u % array(k,iEdge) * 0.5*(state % rho % array(k,cellsOnEdge(1,iEdge)) + state % rho % array(k,cellsOnEdge(2,iEdge)))
+ diag % ru % array(k,iEdge) = state % u % array(k,iEdge) * 0.5*(state % rho_zz % array(k,cellsOnEdge(1,iEdge)) + state % rho_zz % array(k,cellsOnEdge(2,iEdge)))
end do
end do
@@ -3931,11 +3956,11 @@
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
do k=2,grid%nVertLevels
state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
end do
end do
@@ -3947,11 +3972,19 @@
! Calculate surface pressure
do iCell=1,grid%nCells
diag % surface_pressure % array(iCell) = 0.5*gravity/rdzw(1) &
- * (1.25* rho(1,iCell) * (1. + scalars(state % index_qv, 1, iCell)) &
- - 0.25* rho(2,iCell) * (1. + scalars(state % index_qv, 2, iCell)))
+ * (1.25* rho_zz(1,iCell) * (1. + scalars(state % index_qv, 1, iCell)) &
+ - 0.25* rho_zz(2,iCell) * (1. + scalars(state % index_qv, 2, iCell)))
diag % surface_pressure % array(iCell) = diag % surface_pressure % array(iCell) + pp(1,iCell) + ppb(1,iCell)
end do
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(state % index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_gfs
subroutine nhyd_test_case_sst(domain, dminfo, grid, fg, state, diag, test_case, parinfo)
Modified: branches/atmos_physics/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry 2011-08-24 23:29:04 UTC (rev 953)
@@ -162,8 +162,8 @@
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
var persistent real w ( nVertLevelsP1 nCells Time ) 2 iro w state - -
-var persistent real rho ( nVertLevels nCells Time ) 2 iro rho state - -
-var persistent real theta ( nVertLevels nCells Time ) 2 iro theta state - -
+var persistent real rho_zz ( nVertLevels nCells Time ) 2 r rho_zz state - -
+var persistent real theta_m ( nVertLevels nCells Time ) 2 r theta_m state - -
var persistent real qv ( nVertLevels nCells Time ) 2 iro qv state scalars moist
var persistent real qc ( nVertLevels nCells Time ) 2 iro qc state scalars moist
var persistent real qr ( nVertLevels nCells Time ) 2 iro qr state scalars moist
@@ -176,8 +176,8 @@
# Tendency variables
var persistent real tend_u ( nVertLevels nEdges Time ) 1 o u tend - -
var persistent real tend_w ( nVertLevelsP1 nCells Time ) 1 o w tend - -
-var persistent real tend_rho ( nVertLevels nCells Time ) 1 o rho tend - -
-var persistent real tend_theta ( nVertLevels nCells Time ) 1 o theta tend - -
+var persistent real tend_rho ( nVertLevels nCells Time ) 1 o rho_zz tend - -
+var persistent real tend_theta ( nVertLevels nCells Time ) 1 o theta_m tend - -
var persistent real tend_qv ( nVertLevels nCells Time ) 1 o qv tend scalars moist
var persistent real tend_qc ( nVertLevels nCells Time ) 1 o qc tend scalars moist
var persistent real tend_qr ( nVertLevels nCells Time ) 1 o qr tend scalars moist
@@ -200,6 +200,8 @@
var persistent real qv_init ( nVertLevels ) 0 iro qv_init mesh - -
# Diagnostic fields: only written to output
+var persistent real rho ( nVertLevels nCells Time ) 1 io rho diag - -
+var persistent real theta ( nVertLevels nCells Time ) 1 io theta diag - -
var persistent real rh ( nVertLevels nCells Time ) 1 iro rh diag - -
var persistent real v ( nVertLevels nEdges Time ) 1 o v diag - -
var persistent real divergence ( nVertLevels nCells Time ) 1 o divergence diag - -
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-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_mpas_core.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -219,7 +219,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
+ call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % diag, block_ptr % mesh)
block_ptr => block_ptr % next
end do
@@ -229,7 +229,7 @@
end subroutine write_output_frame
- subroutine compute_output_diagnostics(state, grid)
+ subroutine compute_output_diagnostics(state, diag, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields for a domain
!
@@ -244,10 +244,18 @@
implicit none
type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
type (mesh_type), intent(in) :: grid
integer :: i, eoe
- integer :: iEdge, k
+ integer :: iCell, k
+
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * state % scalars % array(state % index_qv,k,iCell))
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * grid % zz % array(k,iCell)
+ end do
+ end do
end subroutine compute_output_diagnostics
Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -135,7 +135,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:), pointer :: surface_pressure
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
- real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
+ real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt
real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3, zb, zb3
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
@@ -243,9 +243,9 @@
ppb => diag % pressure_base % array
pp => diag % pressure_p % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
surface_pressure => diag % surface_pressure % array
@@ -629,7 +629,7 @@
p (k,i) = ((ppb(k,i)+pp(k,i))/p0)**(rgas/cp)
t (k,i) = tt(k)/p(k,i)
rt (k,i) = t(k,i)*rr(k,i)+rb(k,i)*(t(k,i)-tb(k,i))
- rho (k,i) = rb(k,i) + rr(k,i)
+ rho_zz (k,i) = rb(k,i) + rr(k,i)
end do
!calculation of surface pressure:
@@ -698,7 +698,7 @@
cell1 = grid % CellsOnEdge % array(1,iEdge)
cell2 = grid % CellsOnEdge % array(2,iEdge)
do k=1,nz1
- diag % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+ diag % ru % array (k,iEdge) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*state % u % array (k,iEdge)
end do
!
@@ -795,11 +795,11 @@
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
do k=2,grid%nVertLevels
state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
end do
end do
@@ -827,6 +827,14 @@
write(0,*) ' i, psurf, lat ',i,psurf,grid%latCell%array(i)*180./3.1415828
enddo
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_jw
subroutine calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
@@ -989,7 +997,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
!This is temporary variable here. It just need when calculate tangential velocity v.
@@ -1068,11 +1076,11 @@
p => diag % exner % array
cqw => diag % cqw % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
pp => diag % pressure_p % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
u => state % u % array
ru => diag % ru % array
@@ -1439,7 +1447,7 @@
do i=1,grid % ncells
do k=1,nz1
grid % t_init % array(k,i) = t_init_1d(k)
- rho(k,i) = rb(k,i)+rr(k,i)
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
end do
end do
@@ -1448,7 +1456,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- ru (k,i) = 0.5*(rho(k,cell1)+rho(k,cell2))*u(k,i)
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
end do
end if
end do
@@ -1497,6 +1505,14 @@
! write(0,'(i2,3(2x,f14.10)') k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
! end do
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_squall_line
@@ -1519,7 +1535,7 @@
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx, cqw
- real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
+ real (kind=RKIND), dimension(:,:), pointer :: ppb, pb, rho_zz, rb, rr, tb, rtb, p, pp, dss, t, rt, u, ru
real (kind=RKIND), dimension(:,:,:), pointer :: scalars, deriv_two, zf, zf3, zb, zb3
!This is temporary variable here. It just need when calculate tangential velocity v.
@@ -1612,11 +1628,11 @@
p => diag % exner % array
cqw => diag % cqw % array
- rho => state % rho % array
+ rho_zz => state % rho_zz % array
pp => diag % pressure_p % array
rr => diag % rho_p % array
- t => state % theta % array
+ t => state % theta_m % array
rt => diag % rtheta_p % array
u => state % u % array
ru => diag % ru % array
@@ -1959,7 +1975,7 @@
do i=1,grid % ncells
do k=1,nz1
- rho(k,i) = rb(k,i)+rr(k,i)
+ rho_zz(k,i) = rb(k,i)+rr(k,i)
end do
do k=1,nz1
@@ -1972,7 +1988,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- ru (k,i) = 0.5*(rho(k,cell1)+rho(k,cell2))*u(k,i)
+ ru (k,i) = 0.5*(rho_zz(k,cell1)+rho_zz(k,cell2))*u(k,i)
end do
end if
end do
@@ -2063,11 +2079,11 @@
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
do k=2,grid%nVertLevels
state % w % array(k,iCell) = diag % rw % array(k,iCell) &
- / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ / (fzp(k) * state % rho_zz % array(k-1,iCell) + fzm(k) * state % rho_zz % array(k,iCell))
end do
end do
@@ -2099,6 +2115,14 @@
! write(0,*) ' k,u_init, t_init, qv_init ',k,grid % u_init % array(k),grid % t_init% array(k),grid % qv_init % array(k)
! end do
+ ! Compute rho and theta from rho_zz and theta_m
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ diag % rho % array(k,iCell) = state % rho_zz % array(k,iCell) * zz(k,iCell)
+ diag % theta % array(k,iCell) = state % theta_m % array(k,iCell) / (1.0 + 1.61 * scalars(index_qv,k,iCell))
+ end do
+ end do
+
end subroutine nhyd_test_case_mtn_wave
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-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -117,8 +117,8 @@
! WCS-parallel: first three and rtheta_p arise from scalar transport and microphysics update (OK). Others come from where?
-! theta
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(1) % state % theta % array(:,:), &
+! theta_m
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(1) % state % theta_m % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
! scalars
@@ -181,7 +181,7 @@
do while (associated(block))
call physics_addtend( domain % dminfo , block % parinfo % cellsToSend, block % parinfo % cellsToRecv , &
block % mesh , block % tend, block % tend_physics , &
- block % state % time_levs(2) % state % rho % array(:,:), &
+ block % state % time_levs(2) % state % rho_zz % array(:,:), &
block % diag % rho_edge % array(:,:) )
block => block % next
end do
@@ -333,7 +333,7 @@
! WCS-parallel: seems like we already use u and w (and other state variables) as if they were already correctly set in the halo,
! but they are not (at least to the outer edges - see communications below? or are those communications redundent?).
-! Perhaps we should communicate u, w, theta, rho, etc after recover_large_step_variables),
+! Perhaps we should communicate u, w, theta_m, rho_zz, etc after recover_large_step_variables),
! cover it with the scalar communications, and then compute solve_diagnostics. I do not think we need to communicate the stuff we compute
! in compute_solve_diagnostics if we compute it out in the halo (and I think we do - the halos should be large enough).
@@ -346,7 +346,7 @@
if(debug) write(0,*) ' diagnostics complete '
- ! need communications here to fill out u, w, theta, p, and pp, scalars, etc
+ ! need communications here to fill out u, w, theta_m, p, and pp, scalars, etc
! so that they are available for next RK step or the first rk substep of the next timestep
block => domain % blocklist
do while (associated(block))
@@ -479,8 +479,8 @@
s_old % u % array = s_new % u % array
s_old % w % array = s_new % w % array
- s_old % theta % array = s_new % theta % array
- s_old % rho % array = s_new % rho % array
+ s_old % theta_m % array = s_new % theta_m % array
+ s_old % rho_zz % array = s_new % rho_zz % array
s_old % scalars % array = s_new % scalars % array
end subroutine rk_integration_setup
@@ -589,7 +589,7 @@
cofwt => diag % cofwt % array
cofrz => diag % cofrz % array
- t => s % theta % array
+ t => s % theta_m % array
dtseps = .5*dts*(1.+epssm)
rcv = rgas/(cp-rgas)
@@ -727,7 +727,7 @@
type (mesh_type) :: grid
real (kind=RKIND), intent(in) :: dts
- real (kind=RKIND), dimension(:,:), pointer :: rho, theta, ru_p, rw_p, rtheta_pp, &
+ real (kind=RKIND), dimension(:,:), pointer :: rho_zz, theta_m, ru_p, rw_p, rtheta_pp, &
rtheta_pp_old, zz, exner, cqu, ruAvg, &
wwAvg, rho_pp, cofwt, coftz, zx, &
a_tri, alpha_tri, gamma_tri, dss, &
@@ -755,8 +755,8 @@
!--
- rho => s % rho % array
- theta => s % theta % array
+ rho_zz => s % rho_zz % array
+ theta_m => s % theta_m % array
w => s % w % array
rtheta_pp => diag % rtheta_pp % array
@@ -780,8 +780,8 @@
dss => grid % dss % array
tend_ru => tend % u % array
- tend_rho => tend % rho % array
- tend_rt => tend % theta % array
+ tend_rho => tend % rho_zz % array
+ tend_rt => tend % theta_m % array
tend_rw => tend % w % array
zz => grid % zz % array
@@ -863,7 +863,7 @@
rs(k,cell1) = rs(k,cell1)-flux/AreaCell(cell1)
rs(k,cell2) = rs(k,cell2)+flux/AreaCell(cell2)
- flux = flux*0.5*(theta(k,cell2)+theta(k,cell1))
+ 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)
@@ -914,7 +914,7 @@
do k=2,nVertLevels
rw_p(k,iCell) = (rw_p(k,iCell)-dts*dss(k,iCell)* &
(fzm(k)*zz (k,iCell)+fzp(k)*zz (k-1,iCell)) &
- *(fzm(k)*rho(k,iCell)+fzp(k)*rho(k-1,iCell)) &
+ *(fzm(k)*rho_zz(k,iCell)+fzp(k)*rho_zz(k-1,iCell)) &
*w(k,iCell) )/(1.+dts*dss(k,iCell))
wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.+epssm)*rw_p(k,iCell)
@@ -945,9 +945,9 @@
real (kind=RKIND), dimension(:,:), pointer :: wwAvg, rw_save, w, rw, rw_p, rtheta_p, rtheta_pp, &
rtheta_p_save, rt_diabatic_tend, rho_p, rho_p_save, &
- rho_pp, rho, rho_base, ruAvg, ru_save, ru_p, u, ru, &
+ rho_pp, rho_zz, rho_base, ruAvg, ru_save, ru_p, u, ru, &
exner, exner_base, rtheta_base, pressure_p, &
- zz, theta, pressure_b, qvapor
+ zz, theta_m, pressure_b, qvapor
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell, rdzw, surface_pressure
real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -972,10 +972,10 @@
rtheta_pp => diag % rtheta_pp % array
rtheta_base => diag % rtheta_base % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
- theta => s % theta % array
+ theta_m => s % theta_m % array
qvapor => s % scalars % array(s%index_qv,:,:)
- rho => s % rho % array
+ rho_zz => s % rho_zz % array
rho_p => diag % rho_p % array
rho_p_save => diag % rho_p_save % array
rho_pp => diag % rho_pp % array
@@ -1018,7 +1018,7 @@
cf3 = grid % cf3 % scalar
! compute new density everywhere so we can compute u from ru.
- ! we will also need it to compute theta below
+ ! we will also need it to compute theta_m below
do iCell = 1, nCells
@@ -1026,7 +1026,7 @@
if( iCell == 479 ) then
write(0,*) ' k,rho_old,rp_old, rho_pp '
do k=1,nVertLevels
- write(0,*) k, rho(k,iCell) ,rho_p(k,iCell), rho_pp(k,iCell)
+ write(0,*) k, rho_zz(k,iCell) ,rho_p(k,iCell), rho_pp(k,iCell)
enddo
end if
end if
@@ -1035,7 +1035,7 @@
rho_p(k,iCell) = rho_p(k,iCell) + rho_pp(k,iCell)
- rho(k,iCell) = rho_p(k,iCell) + rho_base(k,iCell)
+ rho_zz(k,iCell) = rho_p(k,iCell) + rho_base(k,iCell)
end do
! recover owned-cell values in block
@@ -1062,7 +1062,7 @@
! pick up part of diagnosed w from omega
w(k,iCell) = rw(k,iCell)/( (fzm(k)*zz (k,iCell)+fzp(k)*zz (k-1,iCell)) &
- *(fzm(k)*rho(k,iCell)+fzp(k)*rho(k-1,iCell)) )
+ *(fzm(k)*rho_zz(k,iCell)+fzp(k)*rho_zz(k-1,iCell)) )
end do
w(nVertLevels+1,iCell) = 0.
@@ -1079,13 +1079,13 @@
if (rk_step==3) then
rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) &
- - dt * rho(k,iCell) * rt_diabatic_tend(k,iCell)
+ - dt * rho_zz(k,iCell) * rt_diabatic_tend(k,iCell)
else
rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell)
end if
- theta(k,iCell) = (rtheta_p(k,iCell) + rtheta_base(k,iCell))/rho(k,iCell)
+ theta_m(k,iCell) = (rtheta_p(k,iCell) + rtheta_base(k,iCell))/rho_zz(k,iCell)
exner(k,iCell) = (zz(k,iCell)*(rgas/p0)*(rtheta_p(k,iCell)+rtheta_base(k,iCell)))**rcv
! pressure below is perturbation pressure
pressure_p(k,iCell) = zz(k,iCell) * rgas * (exner(k,iCell)*rtheta_p(k,iCell)+rtheta_base(k,iCell) &
@@ -1096,8 +1096,8 @@
!calculation of the surface pressure:
surface_pressure(iCell) = 0.5*gravity/rdzw(1) &
- * (1.25* rho(1,iCell) * (1. + qvapor(1,iCell)) &
- - 0.25* rho(2,iCell) * (1. + qvapor(2,iCell)))
+ * (1.25* rho_zz(1,iCell) * (1. + qvapor(1,iCell)) &
+ - 0.25* rho_zz(2,iCell) * (1. + qvapor(2,iCell)))
surface_pressure(iCell) = surface_pressure(iCell) + pressure_p(1,iCell) + pressure_b(1,iCell)
@@ -1113,7 +1113,7 @@
cell2 = CellsOnEdge(2,iEdge)
! if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed,
- ! though we can avoid division by zero, e.g., by rho(:,cell2)
+ ! though we can avoid division by zero, e.g., by rho_zz(:,cell2)
do k = 1, nVertLevels
@@ -1123,7 +1123,7 @@
ru(k,iEdge) = ru(k,iEdge) + ru_p(k,iEdge)
- u(k,iEdge) = 2.*ru(k,iEdge)/(rho(k,cell1)+rho(k,cell2))
+ u(k,iEdge) = 2.*ru(k,iEdge)/(rho_zz(k,cell1)+rho_zz(k,cell2))
enddo
@@ -1131,26 +1131,26 @@
!SHP-mtn
flux = cf1*ru(1,iEdge) + cf2*ru(2,iEdge) + cf3*ru(3,iEdge)
- w(1,cell2) = w(1,cell2) - zb(1,2,iEdge)*flux/(cf1*rho(1,cell2)+cf2*rho(2,cell2)+cf3*rho(3,cell2))
- w(1,cell1) = w(1,cell1) + zb(1,1,iEdge)*flux/(cf1*rho(1,cell1)+cf2*rho(2,cell1)+cf3*rho(3,cell1))
+ w(1,cell2) = w(1,cell2) - zb(1,2,iEdge)*flux/(cf1*rho_zz(1,cell2)+cf2*rho_zz(2,cell2)+cf3*rho_zz(3,cell2))
+ w(1,cell1) = w(1,cell1) + zb(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
!3rd order stencil
if (config_theta_adv_order == 3) then
w(1,cell2) = w(1,cell2) - sign(1.,ru(1,iEdge))*config_coef_3rd_order &
- *zb3(1,2,iEdge)*flux/(cf1*rho(1,cell2)+cf2*rho(2,cell2)+cf3*rho(3,cell2))
+ *zb3(1,2,iEdge)*flux/(cf1*rho_zz(1,cell2)+cf2*rho_zz(2,cell2)+cf3*rho_zz(3,cell2))
w(2,cell1) = w(2,cell1) + sign(1.,ru(1,iEdge))*config_coef_3rd_order &
- *zb3(1,1,iEdge)*flux/(cf1*rho(1,cell1)+cf2*rho(2,cell1)+cf3*rho(3,cell1))
+ *zb3(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
end if
do k = 2, nVertLevels
flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
- w(k,cell2) = w(k,cell2) - zb(k,2,iEdge)*flux/(fzm(k)*rho(k,cell2)+fzp(k)*rho(k-1,cell2))
- w(k,cell1) = w(k,cell1) + zb(k,1,iEdge)*flux/(fzm(k)*rho(k,cell1)+fzp(k)*rho(k-1,cell1))
+ w(k,cell2) = w(k,cell2) - zb(k,2,iEdge)*flux/(fzm(k)*rho_zz(k,cell2)+fzp(k)*rho_zz(k-1,cell2))
+ w(k,cell1) = w(k,cell1) + zb(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
!3rd order! stencil
if (config_theta_adv_order ==3) then
w(k,cell2) = w(k,cell2) - sign(1.,ru(k,iEdge))*config_coef_3rd_order &
- *zb3(k,2,iEdge)*flux/(fzm(k)*rho(k,cell2)+fzp(k)*rho(k-1,cell2))
+ *zb3(k,2,iEdge)*flux/(fzm(k)*rho_zz(k,cell2)+fzp(k)*rho_zz(k-1,cell2))
w(k,cell1) = w(k,cell1) + sign(1.,ru(k,iEdge))*config_coef_3rd_order &
- *zb3(k,1,iEdge)*flux/(fzm(k)*rho(k,cell1)+fzp(k)*rho(k-1,cell1))
+ *zb3(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
end if
enddo
@@ -1185,7 +1185,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, kdiff
+ real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho_zz, zgrid, kdiff
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, qv_init
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -1225,8 +1225,8 @@
scalar_tend => tend % scalars % array
!**** h_old => s_old % h % array
!**** h_new => s_new % h % array
- h_old => s_old % rho % array
- h_new => s_new % rho % array
+ h_old => s_old % rho_zz % array
+ h_new => s_new % rho_zz % array
wwAvg => diag % wwAvg % array
areaCell => grid % areaCell % array
@@ -1244,7 +1244,7 @@
h_theta_eddy_visc2 = config_h_theta_eddy_visc2
v_theta_eddy_visc2 = config_v_theta_eddy_visc2
rho_edge => diag % rho_edge % array
- rho => s_new % rho % array
+ rho_zz => s_new % rho_zz % array
qv_init => grid % qv_init % array
zgrid => grid % zgrid % array
@@ -1419,7 +1419,7 @@
zp = 0.5*(z3+z4)
do iScalar=1,s_old % num_scalars
- scalar_tend(iScalar,k,iCell) = scalar_tend(iScalar,k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+ scalar_tend(iScalar,k,iCell) = scalar_tend(iScalar,k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
(scalar_new(iScalar,k+1,iCell)-scalar_new(iScalar,k ,iCell))/(zp-z0) &
-(scalar_new(iScalar,k ,iCell)-scalar_new(iScalar,k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -1437,7 +1437,7 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
- scalar_tend(iScalar,k,iCell) = scalar_tend(iScalar,k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+ scalar_tend(iScalar,k,iCell) = scalar_tend(iScalar,k,iCell) + v_theta_eddy_visc2*prandtl*rho_zz(k,iCell)*(&
(-qv_init(k+1)+qv_init(k))/(zp-z0) &
-(-qv_init(k)+qv_init(k-1))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -1568,8 +1568,8 @@
scalar_tend => tend % scalars % array
!**** h_old => s_old % h % array
!**** h_new => s_new % h % array
- h_old => s_old % rho % array
- h_new => s_new % rho % array
+ h_old => s_old % rho_zz % array
+ h_new => s_new % rho_zz % array
wwAvg => diag % wwAvg % array
areaCell => grid % areaCell % array
@@ -1937,8 +1937,8 @@
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, h_theta_eddy_visc4
real (kind=RKIND) :: u_diffusion
real (kind=RKIND), dimension(:), pointer :: fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, meshScalingDel2, meshScalingDel4
- 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, &
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, zgrid, rho_edge, rho_zz, ru, u, v, tend_u, &
+ circulation, divergence, vorticity, ke, pv_edge, theta_m, rw, tend_rho, &
rt_diabatic_tend, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
h_divergence, kdiff
@@ -1998,7 +1998,7 @@
coef_3rd_order = config_coef_3rd_order
- rho => s % rho % array
+ rho_zz => s % rho_zz % array
rho_edge => diag % rho_edge % array
rb => diag % rho_base % array
rr => diag % rho_p % array
@@ -2008,7 +2008,7 @@
ru => diag % ru % array
w => s % w % array
rw => diag % rw % array
- theta => s % theta % array
+ theta_m => s % theta_m % array
circulation => diag % circulation % array
divergence => diag % divergence % array
vorticity => diag % vorticity % array
@@ -2042,9 +2042,9 @@
tend_u => tend % u % array
- tend_theta => tend % theta % array
+ tend_theta => tend % theta_m % array
tend_w => tend % w % array
- tend_rho => tend % rho % array
+ tend_rho => tend % rho_zz % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
tend_u_euler => tend % u_euler % array
@@ -2740,11 +2740,11 @@
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels
-! tend_w(k,iCell) = tend_w(k,iCell) + v_mom_eddy_visc2*0.5*(rho(k,iCell)+rho(k-1,iCell))*( &
+! tend_w(k,iCell) = tend_w(k,iCell) + v_mom_eddy_visc2*0.5*(rho_zz(k,iCell)+rho_zz(k-1,iCell))*( &
! (w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
! -(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
- tend_w_euler(k,iCell) = tend_w_euler(k,iCell) + v_mom_eddy_visc2*0.5*(rho(k,iCell)+rho(k-1,iCell))*( &
+ tend_w_euler(k,iCell) = tend_w_euler(k,iCell) + v_mom_eddy_visc2*0.5*(rho_zz(k,iCell)+rho_zz(k-1,iCell))*( &
(w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
-(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
end do
@@ -2771,9 +2771,9 @@
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels
tend_w(k,iCell) = tend_w(k,iCell) &
- + rho(k,iCell)*( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
+ + rho_zz(k,iCell)*( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
+(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
- + 2.*omega_e*cos(grid % latCell % array(iCell))*rho(k,iCell) &
+ + 2.*omega_e*cos(grid % latCell % array(iCell))*rho_zz(k,iCell) &
*(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))
end do
@@ -2796,7 +2796,7 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,grid % nVertLevels
- flux = dvEdge(iEdge) * ru(k,iEdge) * ( 0.5*(theta(k,cell1) + theta(k,cell2)) )
+ flux = dvEdge(iEdge) * ru(k,iEdge) * ( 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) )
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
tend_theta(k,cell2) = tend_theta(k,cell2) + flux
end do
@@ -2812,35 +2812,35 @@
do k=1,grid % nVertLevels
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta_m(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta_m(k,cell2)
do i=1, grid % nEdgesOnCell % array (cell1)
if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell1))
end do
do i=1, grid % nEdgesOnCell % array (cell2)
if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell2))
end do
! 3rd order stencil
if( u(k,iEdge) > 0) then
flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
else
flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
! flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
-! 0.5*(theta(k,cell1) + theta(k,cell2)) &
+! 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) &
! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
! else
! flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
-! 0.5*(theta(k,cell1) + theta(k,cell2)) &
+! 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) &
! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
end if
@@ -2860,19 +2860,19 @@
do k=1,grid % nVertLevels
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta_m(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta_m(k,cell2)
do i=1, grid % nEdgesOnCell % array (cell1)
if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
- d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell1))
end do
do i=1, grid % nEdgesOnCell % array (cell2)
if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
- d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell2))
end do
flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
@@ -2885,7 +2885,7 @@
end if
!
- ! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
+ ! horizontal mixing for theta_m - 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)
!
@@ -2902,7 +2902,7 @@
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)
+ theta_turb_flux = h_theta_eddy_visc2*prandtl*(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
@@ -2923,7 +2923,7 @@
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)
+ *(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
@@ -2948,8 +2948,8 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
do k=1,grid % nVertLevels
- delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*rho_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*rho_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*rho_edge(k,iEdge)*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
+ delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*rho_edge(k,iEdge)*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
end do
end do
@@ -2995,28 +2995,28 @@
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))
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta_m(k,iCell)+fzp(k)*theta_m(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))
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta_m(k,iCell)+fzp(k)*theta_m(k-1,iCell))
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 )
+ wdtz(k) = flux3( theta_m(k-2,iCell),theta_m(k-1,iCell),theta_m(k,iCell),theta_m(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))
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta_m(k,iCell)+fzp(k)*theta_m(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))
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta_m(k,iCell)+fzp(k)*theta_m(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) )
+ wdtz(k) = flux4( theta_m(k-2,iCell),theta_m(k-1,iCell),theta_m(k,iCell),theta_m(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))
+ wdtz(k) = rw(k,icell)*(fzm(k)*theta_m(k,iCell)+fzp(k)*theta_m(k-1,iCell))
end if
@@ -3025,7 +3025,7 @@
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,iCell) + rho(k,iCell)*rt_diabatic_tend(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + rho_zz(k,iCell)*rt_diabatic_tend(k,iCell)
end do
end do
@@ -3050,12 +3050,12 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
-! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
-! (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
-! -(theta(k ,iCell)-theta(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(k,iCell)*(&
- (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
- -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*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)*(&
+ (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
end do
@@ -3072,12 +3072,12 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
-! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
-! ((theta(k+1,iCell)-t_init(k+1,iCell))-(theta(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
-! -((theta(k ,iCell)-t_init(k,iCell))-(theta(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(k,iCell)*(&
- ((theta(k+1,iCell)-t_init(k+1,iCell))-(theta(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
- -((theta(k ,iCell)-t_init(k,iCell))-(theta(k-1,iCell)-t_init(k-1,iCell)))/(z0-zm) )/(0.5*(zp-zm))
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*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)*(&
+ ((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
end do
@@ -3131,7 +3131,7 @@
real (kind=RKIND) :: ke_fact
- h => s % rho % array
+ h => s % rho_zz % array
u => s % u % array
v => diag % v % array
vh => diag % rv % array
@@ -3399,27 +3399,34 @@
rcv = rgas / (cp-rgas)
p0 = 1.e5 ! this should come from somewhere else...
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ state % theta_m % array(k,iCell) = diag % theta % array(k,iCell) * (1.0 + 1.61 * state % scalars % array(state % index_qv,k,iCell))
+ state % rho_zz % array(k,iCell) = diag % rho % array(k,iCell) / grid % zz % array(k,iCell)
+ end do
+ end do
+
do iEdge = 1, grid % nEdges
iCell1 = grid % cellsOnEdge % array(1,iEdge)
iCell2 = grid % cellsOnEdge % array(2,iEdge)
do k=1,grid % nVertLevels
- diag % ru % array(k,iEdge) = 0.5 * state % u % array(k,iEdge) * (state % rho % array(k,iCell1) + state % rho % array(k,iCell2))
+ diag % ru % array(k,iEdge) = 0.5 * state % u % array(k,iEdge) * (state % rho_zz % array(k,iCell1) + state % rho_zz % array(k,iCell2))
end do
end do
- ! Compute w from rho and rw
+ ! Compute w from rho_zz and rw
do iCell=1,grid%nCells
diag % rw % array(1,iCell) = 0.
diag % rw % array(grid%nVertLevels+1,iCell) = 0.
do k=2,grid%nVertLevels
diag % rw % array(k,iCell) = state % w % array(k,iCell) &
- * (grid % fzp % array(k) * state % rho % array(k-1,iCell) + grid % fzm % array(k) * state % rho % array(k,iCell))
+ * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell))
end do
end do
do iCell = 1, grid % nCells
do k=1,grid % nVertLevels
- diag % rho_p % array(k,iCell) = state % rho % array(k,iCell) - diag % rho_base % array(k,iCell)
+ diag % rho_p % array(k,iCell) = state % rho_zz % array(k,iCell) - diag % rho_base % array(k,iCell)
end do
end do
@@ -3431,8 +3438,8 @@
do iCell = 1, grid % nCells
do k=1,grid % nVertLevels
- diag % rtheta_p % array(k,iCell) = state % theta % array(k,iCell) * diag % rho_p % array(k,iCell) &
- + diag % rho_base % array(k,iCell) * (state % theta % array(k,iCell) - diag % theta_base % array(k,iCell))
+ diag % rtheta_p % array(k,iCell) = state % theta_m % array(k,iCell) * diag % rho_p % array(k,iCell) &
+ + diag % rho_base % array(k,iCell) * (state % theta_m % array(k,iCell) - diag % theta_base % array(k,iCell))
end do
end do
@@ -3481,10 +3488,10 @@
do k = 1, grid % nVertLevels
- tend % rt_diabatic_tend % array(k,iCell) = state_new % theta % array(k,iCell)
+ tend % rt_diabatic_tend % array(k,iCell) = state_new % theta_m % array(k,iCell)
- t(k) = state_new % theta % array(k,iCell)/(1. + 1.61*state_new % scalars % array(state_new % index_qv,k,iCell))
- rho(k) = grid % zz % array(k,iCell)*state_new % rho % array(k,iCell)
+ t(k) = state_new % theta_m % array(k,iCell)/(1. + 1.61*state_new % scalars % array(state_new % index_qv,k,iCell))
+ rho(k) = grid % zz % array(k,iCell)*state_new % rho_zz % array(k,iCell)
p(k) = diag % exner % array(k,iCell)
qv(k) = max(0.,state_new % scalars % array(state_new % index_qv,k,iCell))
qc(k) = max(0.,state_new % scalars % array(state_new % index_qc,k,iCell))
@@ -3499,10 +3506,10 @@
do k = 1, grid % nVertLevels
- state_new % theta % array(k,iCell) = t(k)*(1.+1.61*qv(k))
- tend % rt_diabatic_tend % array(k,iCell) = state_new % rho % array(k,iCell) * &
- (state_new % theta % array(k,iCell) - tend % rt_diabatic_tend % array(k,iCell))/dt
- diag % rtheta_p % array(k,iCell) = state_new % rho % array(k,iCell) * state_new % theta % array(k,iCell) &
+ state_new % theta_m % array(k,iCell) = t(k)*(1.+1.61*qv(k))
+ tend % rt_diabatic_tend % array(k,iCell) = state_new % rho_zz % array(k,iCell) * &
+ (state_new % theta_m % array(k,iCell) - tend % rt_diabatic_tend % array(k,iCell))/dt
+ diag % rtheta_p % array(k,iCell) = state_new % rho_zz % array(k,iCell) * state_new % theta_m % array(k,iCell) &
- diag % rtheta_base % array(k,iCell)
state_new % scalars % array(state_new % index_qv,k,iCell) = qv(k)
state_new % scalars % array(state_new % index_qc,k,iCell) = qc(k)
Modified: branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -143,7 +143,7 @@
real(kind=RKIND),dimension(:),pointer :: fzm,fzp,rdzw
real(kind=RKIND),dimension(:,:),pointer:: zgrid
real(kind=RKIND),dimension(:,:),pointer:: zz,exner,pressure_b,rtheta_p,rtheta_b
- real(kind=RKIND),dimension(:,:),pointer:: rho,theta,qv,pressure_p,u,v,w
+ real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,qv,pressure_p,u,v,w
real(kind=RKIND),dimension(:,:),pointer:: qvs,rh
integer:: ip,iEdg
@@ -175,8 +175,8 @@
rtheta_p => diag % rtheta_p % array
rtheta_b => diag % rtheta_base % array
- rho => state % rho % array
- theta => state % theta % array
+ rho_zz => state % rho_zz % array
+ theta_m => state % theta_m % array
qv => state % scalars % array(state%index_qv,:,:)
w => state % w % array
@@ -207,10 +207,10 @@
v_p(i,k,j) = v(k,i)
zz_p(i,k,j) = zz(k,i)
- rho_p(i,k,j) = zz(k,i) * rho(k,i)
+ rho_p(i,k,j) = zz(k,i) * rho_zz(k,i)
rho_p(i,k,j) = rho_p(i,k,j)*(1.+qv_p(i,k,j))
- th_p(i,k,j) = theta(k,i) / (1. + R_v/R_d * qv(k,i))
- t_p(i,k,j) = theta(k,i) * exner(k,i) / (1. + R_v/R_d * qv(k,i))
+ th_p(i,k,j) = theta_m(k,i) / (1. + R_v/R_d * qv(k,i))
+ t_p(i,k,j) = theta_m(k,i) * exner(k,i) / (1. + R_v/R_d * qv(k,i))
pi_p(i,k,j) = exner(k,i)
pres_p(i,k,j) = pressure_p(k,i) + pressure_b(k,i)
@@ -350,7 +350,7 @@
real(kind=RKIND),dimension(:,:),pointer:: zgrid
real(kind=RKIND),dimension(:,:),pointer:: zz,exner,pressure_b,rtheta_p,rtheta_b
- real(kind=RKIND),dimension(:,:),pointer:: rho,theta,qv,rh,pressure_p
+ real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,qv,rh,pressure_p
real(kind=RKIND),dimension(:,:),pointer:: rt_diabatic_tend
!---------------------------------------------------------------------------------------------
@@ -368,8 +368,8 @@
rtheta_p => diag % rtheta_p % array
rtheta_b => diag % rtheta_base % array
- rho => state % rho % array
- theta => state % theta % array
+ rho_zz => state % rho_zz % array
+ theta_m => state % theta_m % array
rh => diag % rh % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
@@ -386,8 +386,8 @@
do j = jts, jte
do k = kts, kte
do i = its, ite
- rho_p(i,k,j) = zz(k,i) * rho(k,i)
- th_p(i,k,j) = theta(k,i) / (1. + R_v/R_d * max(0.,qv(k,i)))
+ rho_p(i,k,j) = zz(k,i) * rho_zz(k,i)
+ th_p(i,k,j) = theta_m(k,i) / (1. + R_v/R_d * max(0.,qv(k,i)))
pi_p(i,k,j) = exner(k,i)
pres_p(i,k,j) = pressure_b(k,i) + pressure_p(k,i)
@@ -477,7 +477,7 @@
integer:: icount
integer:: i,k,j
real(kind=RKIND),dimension(:,:),pointer:: zz,exner,exner_b,pressure_b,rtheta_p,rtheta_b
- real(kind=RKIND),dimension(:,:),pointer:: rho,theta,pressure_p
+ real(kind=RKIND),dimension(:,:),pointer:: rho_zz,theta_m,pressure_p
real(kind=RKIND),dimension(:,:),pointer:: rt_diabatic_tend
!---------------------------------------------------------------------------------------------
@@ -494,8 +494,8 @@
rtheta_p => diag % rtheta_p % array
rtheta_b => diag % rtheta_base % array
- rho => state % rho % array
- theta => state % theta % array
+ rho_zz => state % rho_zz % array
+ theta_m => state % theta_m % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
@@ -506,13 +506,13 @@
do i = its, ite
!potential temperature and diabatic forcing:
- rt_diabatic_tend(k,i) = theta(k,i)
- theta(k,i) = th_p(i,k,j) * (1. + R_v/R_d * qv_p(i,k,j))
- rt_diabatic_tend(k,i) = (theta(k,i) - rt_diabatic_tend(k,i)) / dt_dyn
+ rt_diabatic_tend(k,i) = theta_m(k,i)
+ theta_m(k,i) = th_p(i,k,j) * (1. + R_v/R_d * qv_p(i,k,j))
+ rt_diabatic_tend(k,i) = (theta_m(k,i) - rt_diabatic_tend(k,i)) / dt_dyn
! rt_diabatic_tend(k,i) = 0.
!density-weigthed perturbation potential temperature:
- rtheta_p(k,i) = rho(k,i) * theta(k,i) - rtheta_b(k,i)
+ rtheta_p(k,i) = rho_zz(k,i) * theta_m(k,i) - rtheta_b(k,i)
!exner function:
exner(k,i) = (zz(k,i)*(R_d/P0)*(rtheta_p(k,i)+rtheta_b(k,i)))**rcv
Modified: branches/atmos_physics/src/core_physics/module_physics_todynamics.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_physics_todynamics.F 2011-08-22 22:28:31 UTC (rev 952)
+++ branches/atmos_physics/src/core_physics/module_physics_todynamics.F 2011-08-24 23:29:04 UTC (rev 953)
@@ -72,7 +72,7 @@
rthratensw => tend_physics % rthratensw % array
tend_u => tend % u % array
- tend_theta => tend % theta % array
+ tend_theta => tend % theta_m % array
tend_scalars => tend % scalars % array
tend_scalars = 0.
</font>
</pre>