<p><b>duda</b> 2011-01-05 18:35:10 -0700 (Wed, 05 Jan 2011)</p><p>BRANCH COMMIT<br>
<br>
- Merge changes to test case initialization from atmos_physics branch.<br>
- Compute surface_pressure diagnostic, and multiply rt_diabatic_tend by rho when computing<br>
tend_theta in compute_dyn_tend() and when computing rtheta_p in recover_large_step_variables(),<br>
as is done in the atmos_physics branch.<br>
<br>
With these changes, results between atmos_physics and atmos_nonhydrostatic branches match with<br>
physics turned off for test case 2.<br>
<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>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2011-01-06 01:04:50 UTC (rev 677)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2011-01-06 01:35:10 UTC (rev 678)
@@ -205,6 +205,9 @@
var persistent real rho_base ( nVertLevels nCells Time ) 1 iro rho_base diag - -
var persistent real theta_base ( nVertLevels nCells Time ) 1 iro theta_base diag - -
+# For moist initialization
+var persistent real qsat ( nVertLevels nCells Time ) 1 o qsat diag_physics - -
+var persistent real relhum ( nVertLevels nCells Time ) 1 o relhum diag_physics - -
var persistent real ruAvg ( nVertLevels nEdges Time ) 1 - ruAvg diag - -
var persistent real wwAvg ( nVertLevelsP1 nCells Time ) 1 - wwAvg diag - -
@@ -243,3 +246,6 @@
# Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 iro coeffs_reconstruct mesh - -
+var persistent real surface_pressure ( nCells Time ) 1 o surface_pressure diag - -
+var persistent real surface_temperature ( nCells Time ) 1 o surface_temperature diag - -
+
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2011-01-06 01:04:50 UTC (rev 677)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2011-01-06 01:35:10 UTC (rev 678)
@@ -44,7 +44,8 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
+ call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, &
+ block_ptr % diag_physics ,config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
@@ -96,7 +97,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_jw(grid, state, diag, test_case)
+ subroutine nhyd_test_case_jw(grid, state, diag, diag_physics, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -106,6 +107,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
type (diag_type), intent(inout) :: diag
+ type (diag_physics_type), intent(inout) :: diag_physics
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: u0 = 35.0
@@ -115,15 +117,23 @@
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
real (kind=RKIND), parameter :: theta_c = pii/4.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rh_max = 0.4 ! Maximum relative humidity
real (kind=RKIND), parameter :: k_x = 9. ! Normal mode wave number
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 :: zf, zf3, zb, zb3
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+
+!.. initialization of moisture:
+ integer:: index_qv
+! real (kind=RKIND),parameter :: rh_max = 0.40 ! Maximum relative humidity
+ real (kind=RKIND),parameter :: rh_max = 0.70 ! Maximum relative humidity
+ real (kind=RKIND),dimension(:,:), pointer:: qsat,relhum
+ real (kind=RKIND),dimension(:,:,:),pointer:: scalars
+ real (kind=RKIND),dimension(grid%nVertLevels,grid%nCells):: qv
+!.. end initialization of moisture.
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
@@ -141,8 +151,8 @@
real (kind=RKIND) :: r_earth, etavs, ztemp, zd, zt, dz, gam, delt, str
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, qv
- real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp
+ real (kind=RKIND), dimension(grid % nVertLevels) :: temperature
+ real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp, ttemp
integer :: iter
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -224,13 +234,25 @@
t => state % theta % array
rt => diag % rtheta_p % array
+ surface_pressure => diag % surface_pressure % array
+
+!.. initialization of moisture:
scalars => state % scalars % array
+ qsat => diag_physics % qsat % array
+ relhum => diag_physics % relhum % array
+ scalars(:,:,:) = 0.0
+ qsat(:,:) = 0.0
+ relhum(:,:) = 0.0
+ qv(:,:) = 0.0
+!.. end initialization of moisture.
- scalars(:,:,:) = 0.
+ surface_pressure(:) = 0.0
call initialize_advection_rk(grid)
call initialize_deformation_weights(grid)
+ index_qv = state % index_qv
+
xnutr = 0.
zd = 12000.
znut = eta_t
@@ -494,15 +516,20 @@
rr (k,i) = 0.
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
- enddo
- end if
-!
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
+! enddo
+! end if
+
+ 200 format(4i6,8(1x,e15.8))
+ 201 format(3i6,8(1x,e15.8))
+ 202 format(2i6,10(1x,e15.8))
+ 203 format(i6,10(1x,e15.8))
+
! iterations to converge temperature as a function of pressure
!
- do itr = 1,10
+ do itr = 1,30
do k=1,nz1
eta (k) = (ppb(k,i)+pp(k,i))/p0
@@ -515,8 +542,8 @@
end do
phi = grid % latCell % array (i)
do k=1,nz1
- tt(k) = 0.
- tt(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
+! tt(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
+ temperature(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
*sqrt(cos(etav(k)))* &
((-2.*sin(phi)**6 &
*(cos(phi)**2+1./3.)+10./63.) &
@@ -524,13 +551,12 @@
+(1.6*cos(phi)**3 &
*(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
-
- !write(0,*) ' k, tt(k) ',k,tt(k)
ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
ptemp = ppb(k,i) + pp(k,i)
-! qv(k,i) = env_qv( ztemp, tt(k), ptemp, 0 )
- qv(k,i) = 0.
+ !.. conversion from virtual temperature to "modified" temperature:
+ tt(k) = temperature(k)*(1+scalars(index_qv,k,i))
+
end do
! do k=2,nz1
! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,i)))
@@ -568,14 +594,31 @@
rho (k,i) = rb(k,i) + rr(k,i)
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
- enddo
- end if
+ !calculation of surface pressure:
+ surface_pressure(i) = 0.5*dzw(1)*gravity &
+ * (1.25*(rr(1,i) + rb(1,i)) * (1. + scalars(index_qv,1,i)) &
+ - 0.25*(rr(2,i) + rb(2,i)) * (1. + scalars(index_qv,2,i)))
+ surface_pressure(i) = surface_pressure(i) + pp(1,i) + ppb(1,i)
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
+! enddo
+! end if
+
end do ! end loop over cells
+ write(0,*)
+ write(0,*) '--- initialization of water vapor:'
+ do iCell = 1, grid % nCells
+ if(iCell == 1 .or. iCell == grid % nCells) then
+ do k = nz1, 1, -1
+ write(0,202) iCell,k,t(k,iCell),relhum(k,iCell),qsat(k,iCell),scalars(index_qv,k,iCell)
+ enddo
+ write(0,*)
+ endif
+ enddo
+
lat_pert = latitude_pert*pii/180.
lon_pert = longitude_pert*pii/180.
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2011-01-06 01:04:50 UTC (rev 677)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2011-01-06 01:35:10 UTC (rev 678)
@@ -937,8 +937,8 @@
rtheta_p_save, rt_diabatic_tend, rho_p, rho_p_save, &
rho_pp, rho, rho_base, ruAvg, ru_save, ru_p, u, ru, &
exner, exner_base, rtheta_base, pressure_p, &
- zz, theta
- real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell
+ zz, theta, 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
@@ -963,6 +963,7 @@
rtheta_base => diag % rtheta_base % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
theta => s % theta % array
+ qvapor => s % scalars % array(s%index_qv,:,:)
rho => s % rho % array
rho_p => diag % rho_p % array
@@ -980,7 +981,10 @@
exner_base => diag % exner_base % array
pressure_p => diag % pressure_p % array
+ pressure_b => diag % pressure_base % array
+ surface_pressure => diag % surface_pressure % array
+ rdzw => grid % rdzw % array
zz => grid % zz % array
zb => grid % zb % array
zb3 => grid % zb3 % array
@@ -1062,7 +1066,8 @@
do k = 1, nVertLevels
if (rk_step==3) then
- rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) - dt * rt_diabatic_tend(k,iCell)
+ rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) &
+ - dt * rho(k,iCell) * rt_diabatic_tend(k,iCell)
else
rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell)
end if
@@ -1077,6 +1082,13 @@
! end if
+ !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)))
+ surface_pressure(iCell) = surface_pressure(iCell) + pressure_p(1,iCell) + pressure_b(1,iCell)
+
+
end do
! recover time-averaged ruAvg on all edges of owned cells (for upcoming scalar transport).
@@ -2903,7 +2915,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) + rt_diabatic_tend(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + rho(k,iCell)*rt_diabatic_tend(k,iCell)
end do
end do
</font>
</pre>