<p><b>duda</b> 2010-09-24 16:57:36 -0600 (Fri, 24 Sep 2010)</p><p>BRANCH COMMIT<br>
<br>
Changes from Sanghun Park:<br>
- Add z-metric term in w equation<br>
- Add mountain-wave test case<br>
- Split squall-line/super-cell test case into test cases 4 and 5, respectively<br>
<br>
<br>
M src/core_nhyd_atmos/mpas_interface.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 src/core_nhyd_atmos/Makefile<br>
A namelist.input.nhyd_atmos_mtn_wave<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_mtn_wave
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_mtn_wave (rev 0)
+++ branches/atmos_nonhydrostatic/namelist.input.nhyd_atmos_mtn_wave 2010-09-24 22:57:36 UTC (rev 511)
@@ -0,0 +1,35 @@
+&nhyd_model
+ config_test_case = 6
+ config_time_integration = 'SRK3'
+ config_dt = 6.
+ config_ntimesteps = 3000
+ config_output_interval = 100
+ config_number_of_sub_steps = 6
+ config_h_mom_eddy_visc2 = 10.
+ config_h_mom_eddy_visc4 = 0.
+ config_v_mom_eddy_visc2 = 10.
+ config_h_theta_eddy_visc2 = 10.
+ config_h_theta_eddy_visc4 = 0.
+ config_v_theta_eddy_visc2 = 10.
+ config_theta_adv_order = 3
+ config_w_adv_order = 3
+ config_scalar_advection = .false.
+ config_positive_definite = .false.
+ config_monotonic = .false.
+ config_mix_full = .false.
+ config_horiz_mixing = '2d_fixed'
+ config_coef_3rd_order = 0.25
+ config_epssm = 0.2
+/
+
+&io
+ config_input_name = 'grid.nc'
+ config_output_name = 'output.nc'
+ config_restart_name = 'restart.nc'
+/
+
+&restart
+ config_restart_interval = 3000
+ config_do_restart = .false.
+ config_restart_time = 1036800.0
+/
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Makefile
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Makefile 2010-09-24 18:32:31 UTC (rev 510)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Makefile 2010-09-24 22:57:36 UTC (rev 511)
@@ -10,7 +10,7 @@
core_hyd: $(OBJS)
ar -ru libdycore.a $(OBJS)
-module_test_cases.o:
+module_test_cases.o: module_advection.o
module_time_integration.o:
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-09-24 18:32:31 UTC (rev 510)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-09-24 22:57:36 UTC (rev 511)
@@ -128,6 +128,10 @@
var persistent real fzp ( nVertLevels ) 0 iro fzp mesh - -
var persistent real zx ( nVertLevelsP1 nEdges ) 0 iro zx mesh - -
var persistent real zz ( nVertLevelsP1 nCells ) 0 iro zz mesh - -
+var persistent real zf ( nVertLevelsP1 TWO nEdges ) 0 iro zf mesh - -
+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 - -
# coefficients for the vertical tridiagonal solve
# Note: these could be local but...
@@ -172,7 +176,7 @@
# var persistent real pp ( nVertLevelsP1 nCells Time ) 2 - pp state - -
var persistent real u_init ( nVertLevels ) 0 iro u_init mesh - -
-var persistent real t_init ( nVertLevels ) 0 iro t_init mesh - -
+var persistent real t_init ( nVertLevels nCells ) 0 iro t_init mesh - -
var persistent real qv_init ( nVertLevels ) 0 iro qv_init mesh - -
# Diagnostic fields: only written to output
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2010-09-24 18:32:31 UTC (rev 510)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2010-09-24 22:57:36 UTC (rev 511)
@@ -4,6 +4,7 @@
use configure
use constants
use dmpar
+ use advection
contains
@@ -46,9 +47,11 @@
block_ptr => block_ptr % next
end do
- else if (config_test_case == 4 ) then
+ else if ((config_test_case == 4) .or. (config_test_case ==5)) then
write(0,*) ' squall line - super cell test case '
+ if (config_test_case == 4) write(0,*) ' squall line test case'
+ if (config_test_case == 5) write(0,*) ' supercell test case'
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
@@ -61,9 +64,25 @@
block_ptr => block_ptr % next
end do
+ else if (config_test_case == 6 ) then
+
+ write(0,*) ' mountain wave test case '
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ write(0,*) ' calling test case setup '
+ call nhyd_test_case_mtn_wave(block_ptr % mesh, block_ptr % state % time_levs(1) % state, 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)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
else
- write(0,*) ' Only test case 1, 2, 3 and 4 are currently supported for nonhydrostatic core '
+
+ write(0,*) ' Only test case 1, 2, 3, 4, 5 and 6 are currently supported for nonhydrostatic core '
stop
end if
@@ -95,14 +114,17 @@
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 :: zf, zf3, zb, zb3
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
- integer, dimension(:,:), pointer :: edgesOnEdge
+ integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
+ real, dimension(:), pointer :: dvEdge, AreaCell
real, dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
@@ -133,6 +155,7 @@
real (kind=RKIND), dimension(nlat, grid % nVertLevels) :: u_2d, etavs_2d
real (kind=RKIND), dimension(nlat) :: lat_2d
real (kind=RKIND) :: dlat
+ real (kind=RKIND) :: z_edge, z_edge3, d2fdx2_cell1, d2fdx2_cell2
!
! Scale all distances and areas from a unit sphere to one with radius a
@@ -155,6 +178,15 @@
weightsOnEdge => grid % weightsOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
edgesOnEdge => grid % edgesOnEdge % array
+ dvEdge => grid % dvEdge % array
+ AreaCell => grid % AreaCell % array
+ CellsOnEdge => grid % CellsOnEdge % array
+
+ deriv_two => grid % deriv_two % array
+ zf => grid % zf % array
+ zf3 => grid % zf3% array
+ zb => grid % zb % array
+ zb3 => grid % zb3% array
nz1 = grid % nVertLevels
nz = nz1 + 1
@@ -188,6 +220,9 @@
scalars(:,:,:) = 0.
+ call initialize_advection_rk(grid)
+ call initialize_deformation_weights(grid)
+
xnutr = 0.
zd = 12000.
znut = eta_t
@@ -579,6 +614,14 @@
state % u % array(k,iEdge) = fluxk + u_pert
end do
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do k=1,nz1
+ grid % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+ end do
+ end if
+
!
! Generate rotated Coriolis field
!
@@ -599,12 +642,89 @@
!
! CALCULATION OF OMEGA, RW = ZX * RU + ZZ * RW
!
- ! we are assuming w and rw are zero for this initialization
- ! i.e., no terrain
+
!
+ ! pre-calculation z-metric terms in omega eqn.
+ !
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else if (config_theta_adv_order == 3 .or. config_theta_adv_order ==4) then !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+ if (k /= 1) then
+ zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb(k,1,iEdge)
+ zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb(k,2,iEdge)
+ zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb3(k,1,iEdge)
+ zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb3(k,2,iEdge)
+ end if
+
+ end do
+
+ end if
+ end do
+
+ ! for including terrain
grid % rw % array = 0.
state % w % array = 0.
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+ do k = 2, grid%nVertLevels
+ flux = (fzm(k)*grid % ru % array(k,iEdge)+fzp(k)*grid % ru % array(k-1,iEdge))
+ grid % rw % array(k,cell2) = grid % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+ grid % rw % array(k,cell1) = grid % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+
+ if (config_theta_adv_order ==3) then
+ grid % rw % array(k,cell2) = grid % rw % array(k,cell2) &
+ - sign(1.,grid % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ grid % rw % array(k,cell1) = grid % rw % array(k,cell1) &
+ + sign(1.,grid % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+ end if
+
+ end do
+
+
!
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
@@ -697,54 +817,39 @@
type (state_type), intent(inout) :: state
integer, intent(in) :: test_case
- real (kind=RKIND), parameter :: u0 = 35.0
- real (kind=RKIND), parameter :: alpha_grid = 0. ! no grid rotation
- real (kind=RKIND), parameter :: omega_e = 7.29212e-05
- real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
- 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 :: 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 :: scalars
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
- integer :: index_qv
-
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
integer, dimension(:), pointer :: nEdgesOnEdge
integer, dimension(:,:), pointer :: edgesOnEdge
real, dimension(:,:), pointer :: weightsOnEdge
- real (kind=RKIND) :: flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, cell1, cell2, nCellsSolve
+ integer :: index_qv
- real (kind=RKIND) :: ptop, p0, phi
- real (kind=RKIND) :: lon_Edge
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znu, znw, znwc, znwv
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znuc, znuv
- real (kind=RKIND) :: r_earth, etavs, ztemp, zd, zt, dz, gam, str
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, rh, thi
- real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp
- integer :: iter
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh, thi, tbi, cqwb
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: znuc, znuv, bn, divh, dpn
+ real (kind=RKIND) :: r, xnutr
+ real (kind=RKIND) :: ztemp, zd, zt, dz, str
- real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
- real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
- real (kind=RKIND), dimension(grid % nVertLevels ) :: eta, etav, teta, ppi, tt
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: qvb
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: t_init_1d
real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2
- real (kind=RKIND) :: ztr, thetar, ttr, thetas, um, us, zts, pitop, ptopb, rcp, rcv
- real (kind=RKIND) :: radx, radz, zcent, xmid, delt, xloc, rad, temp, pres, yloc, ymid, a_scale
+ real (kind=RKIND) :: ztr, thetar, ttr, thetas, um, us, zts, pitop, pibtop, ptopb, ptop, rcp, rcv, p0
+ real (kind=RKIND) :: radx, radz, zcent, xmid, delt, xloc, rad, yloc, ymid, a_scale
+ real (kind=RKIND) :: pres, temp, es, qvs
- index_qv = state % index_qv
-
!
! Scale all distances
!
@@ -804,14 +909,16 @@
scalars => state % scalars % array
+ index_qv = state % index_qv
+
scalars(:,:,:) = 0.
+ call initialize_advection_rk(grid)
+ call initialize_deformation_weights(grid)
+
xnutr = 0.
zd = 12000.
- znut = eta_t
- etavs = (1.-0.252)*pii/2.
- r_earth = a
p0 = 1.e+05
rcp = rgas/cp
rcv = rgas/(cp-rgas)
@@ -833,7 +940,7 @@
zt = 20000.
dz = zt/float(nz1)
- write(0,*) ' dz = ',dz
+! write(0,*) ' dz = ',dz
write(0,*) ' hx computation complete '
do k=1,nz
@@ -860,7 +967,7 @@
! ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
ah(k) = 1.
- write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
+! write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
end do
do k=1,nz1
dzw (k) = zw(k+1)-zw(k)
@@ -923,15 +1030,6 @@
end do
enddo
- do k=1,nz1
- write(0,*) ' k, zgrid(k,1),hx(k,1) ',k,zgrid(k,1),hx(k,1)
- enddo
-
- do k=1,nz1
- write(0,*) ' k, zx(k,1) ',k,zx(k,1)
- enddo
-
- write(0,*) ' grid metrics setup complete '
!
! convective initialization
!
@@ -940,22 +1038,18 @@
ttr = 213.
thetas = 300.5
- write(0,*) ' rgas, cp, gravity ',rgas,cp, gravity
+! write(0,*) ' rgas, cp, gravity ',rgas,cp, gravity
-! no flow
-! um = 0.
-! us = 0.
-! zts = 5000.
-! supercell parameters
+ if ( config_test_case == 4) then ! squall line parameters
+ um = 12.
+ us = 10.
+ zts = 2500.
+ else if (config_test_case == 5) then !supercell parameters
um = 30.
us = 15.
zts = 5000.
-! squall-line parameters
- um = 12.
- us = 10.
- zts = 2500.
+ end if
-
do i=1,grid % nCells
do k=1,nz1
ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
@@ -968,6 +1062,10 @@
if(t(k,i).lt.thetas) t(k,i) = thetas
end if
tb(k,i) = t(k,i)
+ thi(k,i) = t(k,i)
+ tbi(k,i) = t(k,i)
+ cqw(k,i) = 1.
+ cqwb(k,i) = 1.
end do
end do
@@ -996,72 +1094,127 @@
call dmpar_bcast_reals(dminfo, nz1, grid % u_init % array)
!
-! reference sounding based on dry atmosphere
+! for reference sounding
!
- pitop = 1.-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
+ do itr=1,30
+
+ pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
+ pibtop = 1.-.5*dzw(1)*gravity*(1.+qvb(1))/(cp*tb(1,1)*zz(1,1))
do k=2,nz1
- pitop = pitop-dzu(k)*gravity/(cp*.5*(tb(k,1)+tb(k-1,1)) &
+ pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*.5*(t(k,1)+t(k-1,1)) &
*.5*(zz(k,1)+zz(k-1,1)))
-
- write(0,*) k,pitop,tb(k,1),dzu(k),tb(k,1)
+ pibtop = pibtop-dzu(k)*gravity/(cp*cqwb(k,1)*.5*(tb(k,1)+tb(k-1,1)) &
+ *.5*(zz(k,1)+zz(k-1,1)))
+
+ !write(0,*) k,pitop,tb(k,1),dzu(k),tb(k,1)
end do
- pitop = pitop-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
+ pitop = pitop-.5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
+ pibtop = pibtop-.5*dzw(nz1)*gravity*(1.+qvb(nz1))/(cp*tb(nz1,1)*zz(nz1,1))
call dmpar_bcast_real(dminfo, pitop)
+ call dmpar_bcast_real(dminfo, pibtop)
- ptopb = p0*pitop**(1./rcp)
+ ptopb = p0*pibtop**(1./rcp)
write(6,*) 'ptopb = ',.01*ptopb
-
+
do i=1, grid % nCells
- pb(nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*tb(nz1,i)*zz(nz1,i))
- p (nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*t (nz1,i)*zz(nz1,i))
+ pb(nz1,i) = pibtop+.5*dzw(nz1)*gravity*(1.+qvb(nz1))/(cp*tb(nz1,i)*zz(nz1,i))
+ p (nz1,i) = pitop+.5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,i))/(cp*t (nz1,i)*zz(nz1,i))
do k=nz1-1,1,-1
- pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*.5*(tb(k,i)+tb(k+1,i)) &
+ pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*cqwb(k+1,i)*.5*(tb(k,i)+tb(k+1,i)) &
*.5*(zz(k,i)+zz(k+1,i)))
- p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*.5*(t (k,i)+t (k+1,i)) &
+ p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*cqw(k+1,i)*.5*(t (k,i)+t (k+1,i)) &
*.5*(zz(k,i)+zz(k+1,i)))
end do
do k=1,nz1
rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
rtb(k,i) = rb(k,i)*tb(k,i)
rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
- cqw(k,i) = 1.
end do
end do
+ !
+ ! update water vapor mixing ratio from humidity profile
+ !
+ do i= 1,grid%nCells
+ do k=1,nz1
+ temp = p(k,i)*thi(k,i)
+ pres = p0*p(k,i)**(1./rcp)
+ qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
+ scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ end do
+ end do
+
+ do k=1,nz1
+!*********************************************************************
+! QVB = QV INCLUDES MOISTURE IN REFERENCE STATE
+! qvb(k) = scalars(index_qv,k,1)
+
+! QVB = 0 PRODUCES DRY REFERENCE STATE
+ qvb(k) = 0.
+!*********************************************************************
+ end do
+
+ do i= 1,grid%nCells
+ do k=1,nz1
+ t (k,i) = thi(k,i)*(1.+1.61*scalars(index_qv,k,i))
+ tb(k,i) = tbi(k,i)*(1.+1.61*qvb(k))
+ end do
+ do k=2,nz1
+ cqw (k,i) = 1./(1.+.5*(scalars(index_qv,k,i)+scalars(index_qv,k-1,i)))
+ cqwb(k,i) = 1./(1.+.5*(qvb(k)+qvb(k-1)))
+ end do
+ end do
+
+ end do !end of iteration loop
+
write(0,*) ' base state sounding '
+ write(0,*) ' k, pb, rb, tb, rtb, t, rr, p, qvb'
do k=1,grid%nVertLevels
- write(0,*) ' k, pb,rb,tb,rtb,t,rr,p ', k,pb(k,1),rb(k,1),tb(k,1),rtb(k,1),t(k,1),rr(k,1),p(k,1)
+ write (0,'i2,8(2x,f19.15)') k,pb(k,1),rb(k,1),tb(k,1),rtb(k,1),t(k,1),rr(k,1),p(k,1),qvb(k)
end do
-!-------------------------------------------------------------------
-! ITERATIONS TO CONVERGE MOIST SOUNDING
!
+! potential temperature perturbation
+!
! delt = -10.
! delt = -0.01
delt = 3.
radx = 10000.
radz = 1500.
zcent = 1500.
- xmid = 150000.
- ymid = 50000.*cos(pii/6.)
+ if (config_test_case == 4) then ! squall line prameters
+ call dmpar_max_real(dminfo, maxval(grid % xCell % array(:)), xmid)
+ xmid = xmid * 0.5
+ ymid = 0.0 ! Not used for squall line
+ else if (config_test_case == 5) then ! supercell parameters
+ call dmpar_max_real(dminfo, maxval(grid % xCell % array(:)), xmid)
+ call dmpar_max_real(dminfo, maxval(grid % yCell % array(:)), ymid)
+ xmid = xmid * 0.5
+ ymid = ymid * 0.5
+ end if
+
do i=1, grid % nCells
xloc = grid % xCell % array(i) - xmid
- yloc = grid % yCell % array(i) - ymid
- yloc = 0.
-! xloc = 0.
+ if (config_test_case == 4) then
+ yloc = 0. !squall line setting
+ else if (config_test_case == 5) then
+ yloc = grid % yCell % array(i) - ymid !supercell setting
+ end if
+
do k = 1,nz1
- thi(k,i) = t(k,i)
ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
rad =sqrt((xloc/radx)**2+(yloc/radx)**2+((ztemp-zcent)/radz)**2)
if(rad.lt.1) then
- thi(k,i) = t(k,i) + delt*cos(.5*pii*rad)**2
+ thi(k,i) = thi(k,i) + delt*cos(.5*pii*rad)**2
end if
+ t (k,i) = thi(k,i)*(1.+1.61*scalars(index_qv,k,i))
end do
end do
do itr=1,30
+
pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
do k=2,nz1
pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*.5*(t (k,1)+t (k-1,1)) &
@@ -1069,67 +1222,577 @@
end do
pitop = pitop - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
ptop = p0*pitop**(1./rcp)
- write(0,*) 'ptop = ',.01*ptop
+ write(0,*) 'ptop = ',.01*ptop, .01*ptopb
call dmpar_bcast_real(dminfo, pitop)
- do i = 1, grid % nCells
+ do i = 1, grid % nCells
pp(nz1,i) = ptop-ptopb+.5*dzw(nz1)*gravity* &
(rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
do k=nz1-1,1,-1
- pp(k,i) = pp(k+1,i)+.5*dzu(k+1)*gravity* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
- +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+! pp(k,i) = pp(k+1,i)+.5*dzu(k+1)*gravity* &
+! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
+! +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+ pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity*( &
+ fzm(k+1)*(rb(k+1,i)*(scalars(index_qv,k+1,i)-qvb(k+1)) &
+ +rr(k+1,i)*(1.+scalars(index_qv,k+1,i))) &
+ +fzp(k+1)*(rb(k ,i)*(scalars(index_qv,k ,i)-qvb(k)) &
+ +rr(k ,i)*(1.+scalars(index_qv,k ,i))))
end do
+ if (itr==1.and.i==1) then
do k=1,nz1
+ print *, "pp-check", pp(k,i)
+ end do
+ end if
+ do k=1,nz1
rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
-rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
end do
+
+ end do ! loop over cells
+
+ end do ! iteration loop
+!----------------------------------------------------------------------
!
-! update water vapor mixing ratio from humitidty profile
-!
+ do k=1,nz1
+ grid % qv_init % array(k) = scalars(index_qv,k,1)
+ end do
+
+ t_init_1d(:) = t(:,1)
+ call dmpar_bcast_reals(dminfo, nz1, t_init_1d)
+ call dmpar_bcast_reals(dminfo, nz1, grid % qv_init % array)
+
+ 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)
+ end do
+ end do
+
+ do i=1,grid % nEdges
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k=1,nz1
- temp = p(k,i)*thi(k,i)
- pres = p0*p(k,i)**(1./rcp)
- qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
- scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ ru (k,i) = 0.5*(rho(k,cell1)+rho(k,cell2))*u(k,i)
end do
-
- do k=1,nz1
- t (k,i) = thi(k,i)*(1.+1.61*scalars(index_qv,k,i))
- end do
- do k=2,nz1
- cqw(k,i) = 1./(1.+.5*( scalars(index_qv,k-1,i) &
- +scalars(index_qv,k ,i)))
- end do
- end do ! iteration loop
+ end if
+ end do
- end do ! loop over cells
-!----------------------------------------------------------------------
+
+ !
+ ! we are assuming w and rw are zero for this initialization
+ ! i.e., no terrain
+ !
+ grid % rw % array = 0.
+ state % w % array = 0.
+
+ grid % zf % array = 0.
+ grid % zf3% array = 0.
+ grid % zb % array = 0.
+ grid % zb3% array = 0.
+
+ !
+ ! Generate rotated Coriolis field
+ !
+ do iEdge=1,grid % nEdges
+ grid % fEdge % array(iEdge) = 0.
+ end do
+
+ do iVtx=1,grid % nVertices
+ grid % fVertex % array(iVtx) = 0.
+ end do
+
+ !
+ ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
+ !
+ state % v % array(:,:) = 0.0
+ do iEdge = 1, grid%nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ if (eoe > 0) then
+ do k = 1, grid%nVertLevels
+ state % v % array(k,iEdge) = state % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ end do
+ end if
+ end do
+ end do
+
+ ! write(0,*) ' k,u_init, t_init, qv_init '
+ ! do k=1,grid%nVertLevels
+ ! 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
+
+ end subroutine nhyd_test_case_squall_line
+
+
+!----------------------------------------------------------------------------------------------------------
+
+
+ subroutine nhyd_test_case_mtn_wave(grid, state, test_case)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+ integer, intent(in) :: test_case
+
+ real (kind=RKIND), parameter :: t0=288., hm=250.
+
+ 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 :: scalars, deriv_two, zf, zf3, zb, zb3
+
+ !This is temporary variable here. It just need when calculate tangential velocity v.
+ integer :: eoe, j
+ integer, dimension(:), pointer :: nEdgesOnEdge
+ integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge
+ real, dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
+ real, dimension(:,:), pointer :: weightsOnEdge
+
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
+ integer :: index_qv
+
+ real (kind=RKIND) :: ptop, pitop, ptopb, p0, flux, d2fdx2_cell1, d2fdx2_cell2
+
+ real (kind=RKIND) :: ztemp, zd, zt, dz, str
+
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh
+ real (kind=RKIND) :: ptmp, es, qvs, xnutr, ptemp
+ integer :: iter
+
+ real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
+ real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
+
+ real (kind=RKIND) :: d1, d2, d3, cof1, cof2, cf1, cf2, cf3
+ real (kind=RKIND) :: um, us, rcp, rcv
+ real (kind=RKIND) :: xmid, temp, pres, a_scale
+
+ real (kind=RKIND) :: xi, xa, xc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, z_edge, z_edge3
+
+ integer, dimension(grid % nCells, 2) :: next_cell
+ real (kind=RKIND), dimension(grid % nCells) :: hxzt
+ logical, parameter :: terrain_smooth = .false.
+
+ !
+ ! Scale all distances
+ !
+
+ a_scale = 1.0
+
+ grid % xCell % array = grid % xCell % array * a_scale
+ grid % yCell % array = grid % yCell % array * a_scale
+ grid % zCell % array = grid % zCell % array * a_scale
+ grid % xVertex % array = grid % xVertex % array * a_scale
+ grid % yVertex % array = grid % yVertex % array * a_scale
+ grid % zVertex % array = grid % zVertex % array * a_scale
+ grid % xEdge % array = grid % xEdge % array * a_scale
+ grid % yEdge % array = grid % yEdge % array * a_scale
+ grid % zEdge % array = grid % zEdge % array * a_scale
+ grid % dvEdge % array = grid % dvEdge % array * a_scale
+ grid % dcEdge % array = grid % dcEdge % array * a_scale
+ grid % areaCell % array = grid % areaCell % array * a_scale**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * a_scale**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a_scale**2.0
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ dvEdge => grid % dvEdge % array
+ AreaCell => grid % AreaCell % array
+ CellsOnEdge => grid % CellsOnEdge % array
+ deriv_two => grid % deriv_two % array
+
+ nz1 = grid % nVertLevels
+ nz = nz1 + 1
+ nCellsSolve = grid % nCellsSolve
+
+ zgrid => grid % zgrid % array
+ zf => grid % zf % array
+ zf3 => grid % zf3 % array
+ zb => grid % zb % array
+ zb3 => grid % zb3 % array
+ rdzw => grid % rdzw % array
+ dzu => grid % dzu % array
+ rdzu => grid % rdzu % array
+ fzm => grid % fzm % array
+ fzp => grid % fzp % array
+ zx => grid % zx % array
+ zz => grid % zz % array
+ hx => grid % hx % array
+ dss => grid % dss % array
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+
+ ppb => grid % pressure_base % array
+ pb => grid % exner_base % array
+ rb => grid % rho_base % array
+ tb => grid % theta_base % array
+ rtb => grid % rtheta_base % array
+ p => grid % exner % array
+ cqw => grid % cqw % array
+
+ rho => state % rho % array
+
+ pp => state % pressure % array
+ rr => state % rho_p % array
+ t => state % theta % array
+ rt => grid % rtheta_p % array
+ u => state % u % array
+ ru => grid % ru % array
+
+ scalars => state % scalars % array
+
+ index_qv = state % index_qv
+
+ scalars(:,:,:) = 0.
+
+ call initialize_advection_rk(grid)
+ call initialize_deformation_weights(grid)
+
+ xnutr = 0.1
+ zd = 10500.
+
+ p0 = 1.e+05
+ rcp = rgas/cp
+ rcv = rgas/(cp-rgas)
+
+ ! for hx computation
+ xa = 5000. !SHP - should be changed based on grid distance
+ xla = 4000.
+ xc = maxval (grid % xCell % array(:))/2.
+
+ ! metrics for hybrid coordinate and vertical stretching
+ str = 1.0
+ zt = 21000.
+ dz = zt/float(nz1)
+! write(0,*) ' dz = ',dz
+
+ do k=1,nz
+
+! sh(k) is the stretching specified for height surfaces
+
+ zc(k) = zt*(real(k-1)*dz/zt)**str
+
+! to specify specific heights zc(k) for coordinate surfaces,
+! input zc(k)
+! zw(k) is the hieght of zeta surfaces
+! zw(k) = (k-1)*dz yields constant dzeta
+! and nonconstant dzeta/dz
+! zw(k) = sh(k)*zt yields nonconstant dzeta
+! and nearly constant dzeta/dz
+
+! zw(k) = float(k-1)*dz
+ zw(k) = zc(k)
!
- write(0,*) ' sounding for the simulation '
+! ah(k) governs the transition between terrain-following
+! and pureheight coordinates
+! ah(k) = 0 is a terrain-following coordinate
+! ah(k) = 1 is a height coordinate
+
+! ah(k) = 1.-cos(.5*pii*(k-1)*dz/zt)**6
+ ah(k) = 1.
+! write(0,*) ' k, zc, zw, ah ',k,zc(k),zw(k),ah(k)
+ end do
do k=1,nz1
- write(6,10) .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
- .01*p0*p(k,1)**(1./rcp),t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
- 1000.*scalars(index_qv,k,1),u(k,1)
- 10 format(1x,5f10.3)
+ dzw (k) = zw(k+1)-zw(k)
+ rdzw(k) = 1./dzw(k)
+ zu(k ) = .5*(zw(k)+zw(k+1))
+ end do
+ do k=2,nz1
+ dzu (k) = .5*(dzw(k)+dzw(k-1))
+ rdzu(k) = 1./dzu(k)
+ fzp (k) = .5* dzw(k )/dzu(k)
+ fzm (k) = .5* dzw(k-1)/dzu(k)
+ rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
+ rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
+ end do
- grid % t_init % array(k) = t(k,1)
- grid % qv_init % array(k) = scalars(index_qv,k,1)
+!********** how are we storing cf1, cf2 and cf3?
+ d1 = .5*dzw(1)
+ d2 = dzw(1)+.5*dzw(2)
+ d3 = dzw(1)+dzw(2)+.5*dzw(3)
+ !cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ !cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ !cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ cof1 = (2.*dzu(2)+dzu(3))/(dzu(2)+dzu(3))*dzw(1)/dzu(2)
+ cof2 = dzu(2) /(dzu(2)+dzu(3))*dzw(1)/dzu(3)
+ cf1 = fzp(2) + cof1
+ cf2 = fzm(2) - cof1 - cof2
+ cf3 = cof2
+
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
+! setting for terrain
+ do iCell=1,grid % nCells
+ xi = grid % xCell % array(iCell)
+ !====1. for pure cosine mountain
+ ! if(abs(xi-xc).ge.2.*xa) then
+ ! hx(1,iCell) = 0.
+ ! else
+ ! hx(1,iCell) = hm*cos(.5*pii*(xi-xc)/(2.*xa))**2.
+ ! end if
+
+ !====2. for cosine mountain
+ !if(abs(xi-xc).lt.xa) THEN
+ ! hx(1,iCell) = hm*cos(pii*(xi-xc)/xla)**2. *cos(.5*pii*(xi-xc)/xa )**2.
+ ! else
+ ! hx(1,iCell) = 0.
+ ! end if
+
+ !====3. for shock mountain
+ hx(1,iCell) = hm*exp(-((xi-xc)/xa)**2)*cos(pii*(xi-xc)/xla)**2.
+
+ hx(nz,iCell) = zt
+
+!***** SHP -> get the temporary point information for the neighbor cell ->> should be changed!!!!!
+ do i=1,grid % nCells
+ !option 1
+ !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)-sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,1) = i
+ !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)+sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,2) = i
+ !option 2
+ next_cell(iCell,1) = iCell - 8 ! note ny=4
+ next_cell(iCell,2) = iCell + 8 ! note ny=4
+
+ if (xCell(iCell).le. 3.*grid % dcEdge % array(1)) then
+ next_cell(iCell,1) = 1
+ else if (xCell(iCell).ge. maxval(xCell(:))-3.*grid % dcEdge % array(1)) then
+ next_cell(iCell,2) = 1
+ end if
+
+ end do
+ enddo
+
+ write(0,*) ' hx computation complete '
+
+
+! smoothing grid for the upper level >> but not propoer for parallel programing
+ dzmin=.7
+ do k=2,nz1
+ sm = .25*min((zc(k)-zc(k-1))/dz,1.)
+ do i=1,grid % nCells
+ hx(k,i) = hx(k-1,i)
+ end do
+
+ do iter = 1,20 !iteration for smoothing
+
+ do i=1,grid % nCells
+ hxzt(i) = hx(k,i) + sm*(hx(k,next_cell(i,2))-2.*hx(k,i)+hx(k,next_cell(i,1)))
+ end do
+ dzh = zc(k) - zc(k-1)
+ do i=1,grid % nCells
+ dzht = zc(k)+hxzt(i) - zc(k-1)-hx(k-1,i)
+ if(dzht.lt.dzh) dzh = dzht
+ end do
+
+ if(dzh.gt.dzmin*(zc(k)-zc(k-1))) then
+ do i=1,grid % nCells
+ hx(k,i) = hxzt(i)
+ end do
+ else
+ goto 99 !SHP - this algorithm should be changed
+ end if
+
+ end do !end of iteration for smoothing
+99 print *,"PASS-SHP"
end do
- call dmpar_bcast_reals(dminfo, nz1, grid % t_init % array)
- call dmpar_bcast_reals(dminfo, nz1, grid % qv_init % array)
-
+ do iCell=1,grid % nCells
+ do k=1,nz
+ if (terrain_smooth) then
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
+ + (1.-ah(k)) * zc(k)
+ else
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(1,iCell)/zt)+hx(1,iCell)) &
+ + (1.-ah(k)) * zc(k)
+ end if
+ end do
+ do k=1,nz1
+ zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
+ end do
+ end do
+
+ do i=1, grid % nEdges
+ iCell1 = grid % CellsOnEdge % array(1,i)
+ iCell2 = grid % CellsOnEdge % array(2,i)
+ do k=1,nz
+ zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
+ end do
+ end do
+ do i=1, grid % nCells
+ do k=1,nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ dss(k,i) = 0.
+ ztemp = zgrid(k,i)
+ if(ztemp.gt.zd+.1) then
+ dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
+ end if
+ end do
+ enddo
+
+ write(0,*) ' grid metrics setup complete '
+
!
+! mountain wave initialization
+!
+ !SHP-original
+ !zinv = 1000.
+ !SHP-schar case
+ zinv = 3000.
+
+ xn2 = 0.0001
+ xn2m = 0.0000
+ xn2l = 0.0001
+
+ um = 10.
+ us = 0.
+
+ do i=1,grid % nCells
+ do k=1,nz1
+ ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
+ tb(k,i) = t0*(1. + xn2m/gravity*ztemp)
+ if(ztemp .le. zinv) then
+ t (k,i) = t0*(1.+xn2l/gravity*ztemp)
+ else
+ t (k,i) = t0*(1.+xn2l/gravity*zinv+xn2/gravity*(ztemp-zinv))
+ end if
+ rh(k,i) = 0.
+ end do
+ end do
+
+! set the velocity field - we are on a plane here.
+
+ do i=1, grid % nEdges
+ cell1 = grid % CellsOnEdge % array(1,i)
+ cell2 = grid % CellsOnEdge % array(2,i)
+ if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
+ do k=1,nz1
+ ztemp = .25*( zgrid(k,cell1 )+zgrid(k+1,cell1 ) &
+ +zgrid(k,cell2)+zgrid(k+1,cell2))
+ u(k,i) = um
+ if(i == 1 ) grid % u_init % array(k) = u(k,i) - us
+ u(k,i) = cos(grid % angleEdge % array(i)) * (u(k,i) - us)
+ end do
+ end if
+ end do
+
+!
+! reference sounding based on dry atmosphere
+!
+ pitop = 1.-.5*dzw(1)*gravity/(cp*tb(1,1)*zz(1,1))
+ do k=2,nz1
+ pitop = pitop-dzu(k)*gravity/(cp*(fzm(k)*tb(k,1)+fzp(k)*tb(k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop = pitop-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
+ ptopb = p0*pitop**(1./rcp)
+
+ do i=1, grid % nCells
+ pb(nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*tb(nz1,i)*zz(nz1,i))
+ p (nz1,i) = pitop+.5*dzw(nz1)*gravity/(cp*t (nz1,i)*zz(nz1,i))
+ do k=nz1-1,1,-1
+ pb(k,i) = pb(k+1,i) + dzu(k+1)*gravity/(cp*.5*(tb(k,i)+tb(k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ p (k,i) = p (k+1,i) + dzu(k+1)*gravity/(cp*.5*(t (k,i)+t (k+1,i)) &
+ *.5*(zz(k,i)+zz(k+1,i)))
+ end do
+ do k=1,nz1
+ rb (k,i) = pb(k,i)**(1./rcv)/((rgas/p0)*tb(k,i)*zz(k,i))
+ rtb(k,i) = rb(k,i)*tb(k,i)
+ rr (k,i) = p (k,i)**(1./rcv)/((rgas/p0)*t (k,i)*zz(k,i))-rb(k,i)
+ cqw(k,i) = 1.
+ end do
+ end do
+
+ write(0,*) ' ***** base state sounding ***** '
+ write(0,*) 'k pb p rb rtb rr tb t'
+ do k=1,grid%nVertLevels
+ write(0,'i2,7(2x,f14.9)') k,pb(k,1),p(k,1),rb(k,1),rtb(k,1),rr(k,1),tb(k,1),t(k,1)
+ end do
+
+ scalars(index_qv,:,:) = 0.
+
+!-------------------------------------------------------------------
+! ITERATIONS TO CONVERGE MOIST SOUNDING
+ do itr=1,30
+ pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1))
+
+ do k=2,nz1
+ pitop = pitop-dzu(k)*gravity/(cp*cqw(k,1)*(fzm(k)*t (k,1)+fzp(k)*t (k-1,1)) &
+ *(fzm(k)*zz(k,1)+fzp(k)*zz(k-1,1)))
+ end do
+ pitop = pitop - .5*dzw(nz1)*gravity*(1.+scalars(index_qv,nz1,1))/(cp*t(nz1,1)*zz(nz1,1))
+ ptop = p0*pitop**(1./rcp)
+
+ do i = 1, grid % nCells
+
+ pp(nz1,i) = ptop-ptopb+.5*dzw(nz1)*gravity* &
+ (rr(nz1,i)+(rr(nz1,i)+rb(nz1,i))*scalars(index_qv,nz1,i))
+ do k=nz1-1,1,-1
+ pp(k,i) = pp(k+1,i)+dzu(k+1)*gravity* &
+ (fzm(k)*(rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i)) &
+ +fzp(k)*(rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i)))
+ end do
+ do k=1,nz1
+ rt(k,i) = (pp(k,i)/(rgas*zz(k,i)) &
+ -rtb(k,i)*(p(k,i)-pb(k,i)))/p(k,i)
+ p (k,i) = (zz(k,i)*(rgas/p0)*(rtb(k,i)+rt(k,i)))**rcv
+ rr(k,i) = (rt(k,i)-rb(k,i)*(t(k,i)-tb(k,i)))/t(k,i)
+ end do
+!
+! update water vapor mixing ratio from humitidty profile
+!
+ do k=1,nz1
+ temp = p(k,i)*t(k,i)
+ pres = p0*p(k,i)**(1./rcp)
+ qvs = 380.*exp(17.27*(temp-273.)/(temp-36.))/pres
+ scalars(index_qv,k,i) = amin1(0.014,rh(k,i)*qvs)
+ end do
+
+ do k=1,nz1
+ t (k,i) = t(k,i)*(1.+1.61*scalars(index_qv,k,i))
+ end do
+ do k=2,nz1
+ cqw(k,i) = 1./(1.+.5*( scalars(index_qv,k-1,i) &
+ +scalars(index_qv,k ,i)))
+ end do
+
+ end do ! loop over cells
+
+ end do ! iteration loop
+!----------------------------------------------------------------------
+!
+ write(0,*) ' *** sounding for the simulation ***'
+ write(0,*) ' z theta pres qv rho_m u rr'
+ do k=1,nz1
+ write(6,'8(f14.9,2x)') .5*(zgrid(k,1)+zgrid(k+1,1))/1000., &
+ t(k,1)/(1.+1.61*scalars(index_qv,k,1)), &
+ .01*p0*p(k,1)**(1./rcp), &
+ 1000.*scalars(index_qv,k,1), &
+ (rb(k,1)+rr(k,1))*(1.+scalars(index_qv,k,1)), &
+ grid % u_init % array(k), rr(k,1)
+ end do
+
do i=1,grid % ncells
do k=1,nz1
rho(k,i) = rb(k,i)+rr(k,i)
end do
+
+ do k=1,nz1
+ grid % t_init % array(k,i) = t(k,i)
+ end do
end do
do i=1,grid % nEdges
@@ -1143,33 +1806,92 @@
end do
!
-! CALCULATION OF OMEGA, RW = ZX * RU + ZZ * RW
+! pre-calculation z-metric terms in omega eqn.
!
-! we are assuming w and rw are zero for this initialization
-! i.e., no terrain
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+ if (k /= 1) then
+ zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
+ zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
+ zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
+ zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
+ end if
+
+ end do
+
+ end if
+ end do
+
+! for including terrain
+ state % w % array(:,:) = 0.0
+ grid % rw % array(:,:) = 0.0
+
!
- grid % rw % array = 0.
- state % w % array = 0.
+! calculation of omega, rw = zx * ru + zz * rw
+!
-! DO I=1,NX
-! IM1=I-1
-! IF(IPER.EQ.1.AND.I.EQ.1) IM1=NX1
-! RW(1 ,I) = 0.
-! RW(NZ,I) = 0.
-! DO K=2,NZ1
-! RW(K ,I) = (FZM(K)*ZZ(K,I)+FZP(K)*ZZ(K-1,I))*(
-! & -RDX*(RUZ(K,I )*(ZUW(K,I )-ZGRID(K,I))
-! & -RUZ(K,IM1)*(ZUW(K,IM1)-ZGRID(K,I))))
-! END DO
-! DO K=1,NZ
-! RW1(K,I) = RW(K,I)
-! END DO
-! END DO
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
- !
- ! Generate rotated Coriolis field
- !
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+ do k = 2, grid%nVertLevels
+ flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
+ grid % rw % array(k,cell2) = grid % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+ grid % rw % array(k,cell1) = grid % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+
+ if (config_theta_adv_order ==3) then
+ grid % rw % array(k,cell2) = grid % rw % array(k,cell2) &
+ - sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ grid % rw % array(k,cell1) = grid % rw % array(k,cell1) &
+ + sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+ end if
+
+ end do
+
+
do iEdge=1,grid % nEdges
grid % fEdge % array(iEdge) = 0.
end do
@@ -1193,16 +1915,14 @@
end do
end do
-! do iCell = 1, grid % nCells
-! rt(5,iCell) = rt(5,iCell) + .1
-! enddo
+! do k=1,grid%nVertLevels
+! 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
+ end subroutine nhyd_test_case_mtn_wave
- do k=1,grid%nVertLevels
- 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
- end subroutine nhyd_test_case_squall_line
+!----------------------------------------------------------------------------------------------------------
real function sphere_distance(lat1, lon1, lat2, lon2, radius)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-09-24 18:32:31 UTC (rev 510)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-09-24 22:57:36 UTC (rev 511)
@@ -239,9 +239,10 @@
block => domain % blocklist
do while (associated(block))
- call recover_large_step_variables( block % state % time_levs(2) % state, &
- block % mesh, rk_sub_timestep(rk_step), &
- number_sub_steps(rk_step) )
+ call recover_large_step_variables( block % state % time_levs(2) % state, &
+ block % mesh, rk_timestep(rk_step), &
+ number_sub_steps(rk_step), rk_step )
+
block => block % next
end do
@@ -626,8 +627,20 @@
type (state_type) :: s_new, s_old
type (tend_type) :: tend
type (mesh_type) :: grid
- integer :: iCell, k
+ integer :: iCell, iEdge, k, cell1, cell2
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ real, dimension(:,:,:), pointer :: zf, zf3
+ real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
+ real (kind=RKIND) :: flux
+ zf => grid % zf % array
+ zf3 => grid % zf3 % array
+ fzm => grid % fzm % array
+ fzp => grid % fzp % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+
grid % rho_pp % array = grid % rho_p_save % array - s_new % rho_p % array
grid % ru_p % array = grid % ru_save % array - grid % ru % array
@@ -637,12 +650,33 @@
do iCell = 1, grid % nCellsSolve
do k = 2, grid % nVertLevels
- tend % w % array(k,iCell) = ( grid % fzm % array (k) * grid % zz % array(k ,iCell) + &
- grid % fzp % array (k) * grid % zz % array(k-1,iCell) ) &
+ tend % w % array(k,iCell) = ( fzm(k) * grid % zz % array(k ,iCell) + &
+ fzp(k) * grid % zz % array(k-1,iCell) ) &
* tend % w % array(k,iCell)
end do
end do
+ do iEdge = 1,grid % nEdges
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k = 2, grid%nVertLevels
+ flux = fzm(k) * tend % u % array(k,iEdge) + fzp(k) * tend % u % array(k-1,iEdge)
+ tend % w % array(k,cell2) = tend % w % array(k,cell2) + zf(k,2,iEdge)*flux
+ tend % w % array(k,cell1) = tend % w % array(k,cell1) - zf(k,1,iEdge)*flux
+!3rd order stencil
+ if (config_theta_adv_order == 3) then
+ tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.,tend % u % array(k,iEdge)) &
+ *config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.,tend % u % array(k,iEdge)) &
+ *config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+
+ end do
+
grid % ruAvg % array = 0.
grid % wwAvg % array = 0.
@@ -671,7 +705,6 @@
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: dpzx
real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells+1 ) :: ts, rs
- real (kind=RKIND), dimension( grid % nVertLevels + 1 , grid % nCells+1 ) :: ws
integer :: cell1, cell2, iEdge, iCell, k
real (kind=RKIND) :: pgrad, flux1, flux2, flux, resm, epssm
@@ -747,7 +780,6 @@
ts = 0.
rs = 0.
- ws = 0.
! acoustic step divergence damping - forward weight rtheta_pp
rtheta_pp_old = rtheta_pp + smdiv*(rtheta_pp - rtheta_pp_old)
@@ -805,14 +837,6 @@
end do
- do k=2,nVertLevels
- flux = 0.5*dvEdge(iEdge)*((zgrid(k,cell2)-zgrid(k,cell1))*(fzm(k)*du(k)+fzp(k)*du(k-1)) )
- flux2 = - (fzm(k)*zz(k ,cell2) +fzp(k)*zz(k-1,cell2))*flux/AreaCell(cell2)
- flux1 = - (fzm(k)*zz(k ,cell1) +fzp(k)*zz(k-1,cell1))*flux/AreaCell(cell1)
- ws(k,cell2) = ws(k,cell2) + flux2
- ws(k,cell1) = ws(k,cell1) + flux1
- enddo
-
end if ! end test for block-owned cells
end do ! end loop over edges
@@ -834,7 +858,7 @@
wwavg(k,iCell) = wwavg(k,iCell) + 0.5*(1.-epssm)*rw_p(k,iCell)
- rw_p(k,iCell) = rw_p(k,iCell) + ws(k,iCell) + dts*tend_rw(k,iCell) &
+ rw_p(k,iCell) = rw_p(k,iCell) + dts*tend_rw(k,iCell) &
- cofwz(k,iCell)*((zz(k ,iCell)*ts (k ,iCell) &
-zz(k-1,iCell)*ts (k-1,iCell)) &
+resm*(zz(k ,iCell)*rtheta_pp(k ,iCell) &
@@ -875,21 +899,22 @@
!------------------------
- subroutine recover_large_step_variables( s, grid, dt, ns )
+ subroutine recover_large_step_variables( s, grid, dt, ns, rk_step )
implicit none
type (state_type) :: s
type (mesh_type) :: grid
- integer, intent(in) :: ns
+ integer, intent(in) :: ns, rk_step
real (kind=RKIND), intent(in) :: dt
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, &
exner, exner_base, rtheta_base, pressure_p, &
- zz, theta, zgrid
+ zz, theta
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell
- integer, dimension(:,:), pointer :: CellsOnEdge
+ real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
+ integer, dimension(:,:), pointer :: cellsOnEdge
integer :: iCell, iEdge, k, cell1, cell2
integer :: nVertLevels, nCells, nCellsSolve, nEdges, nEdgesSolve
@@ -931,7 +956,8 @@
pressure_p => s % pressure % array
zz => grid % zz % array
- zgrid => grid % zgrid % array
+ zb => grid % zb % array
+ zb3 => grid % zb3 % array
fzm => grid % fzm % array
fzp => grid % fzp % array
dvEdge => grid % dvEdge % array
@@ -1009,7 +1035,11 @@
do k = 1, nVertLevels
- rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) ! - dt * rt_diabatic_tend(k,iCell)
+ if (rk_step==3) then
+ rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) - dt * 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)
@@ -1042,14 +1072,29 @@
u(k,iEdge) = 2.*ru(k,iEdge)/(rho(k,cell1)+rho(k,cell2))
enddo
- flux = dvEdge(iEdge)*0.5*(cf1*u(1,iEdge)+cf2*u(2,iEdge)+cf3*u(3,iEdge))*(zgrid(1,cell2)-zgrid(1,cell1))
- w(1,cell2) = w(1,cell2)+flux/AreaCell(cell2)
- w(1,cell1) = w(1,cell1)+flux/AreaCell(cell1)
+ !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))
+!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))
+ 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))
+ end if
do k = 2, nVertLevels
- flux = dvEdge(iEdge)*0.5*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))*(zgrid(k,cell2)-zgrid(k,cell1))
- w(k,cell2) = w(k,cell2)+flux/AreaCell(cell2)
- w(k,cell1) = w(k,cell1)+flux/AreaCell(cell1)
+ 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))
+!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))
+ 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))
+ end if
enddo
! end if
@@ -1812,7 +1857,8 @@
real (kind=RKIND) :: theta_edge, theta_turb_flux, z1, z2, z3, z4, zm, z0, zp, r
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, pgrad
- real (kind=RKIND), dimension(:), pointer :: rdzu, rdzw, fzm, fzp, t_init
+ real (kind=RKIND), dimension(:), pointer :: rdzu, rdzw, fzm, fzp, qv_init
+ real (kind=RKIND), dimension(:,:), pointer :: t_init
real (kind=RKIND), allocatable, dimension(:,:) :: rv, divergence_ru, qtot
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_theta, delsq_divergence
@@ -1892,6 +1938,7 @@
h_diabatic => grid % rt_diabatic_tend % array
t_init => grid % t_init % array
+ qv_init => grid % qv_init % array
rdzu => grid % rdzu % array
rdzw => grid % rdzw % array
@@ -2104,12 +2151,6 @@
wduz(nVertLevels+1) = 0.
- wduz(1) = 0.
- do k=2,nVertLevels
- wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
- end do
- wduz(nVertLevels+1) = 0.
-
do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) - rdzw(k)*(wduz(k+1)-wduz(k))
end do
@@ -2553,23 +2594,25 @@
do k=2,nVertLevels
tend_w(k,iCell) = tend_w(k,iCell)/areaCell(iCell) -rdzu(k)*(wdwz(k+1)-wdwz(k)) &
+!SHP-w
- cqw(k,iCell)*( rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+ gravity* &
-!shpark
- ( fzm(k)*rr(k,iCell) + fzm(k)*(rb(k,iCell)+rr(k,iCell))*qtot(k,iCell) &
- +fzp(k)*rr(k-1,iCell) + fzp(k)*(rb(k-1,iCell)+rr(k-1,iCell))*qtot(k-1,iCell) ))
-
+ ( fzm(k)*(rb(k,iCell)*(qtot(k,iCell)) + &
+ rr(k,iCell)*(1.+qtot(k,iCell))) &
+ +fzp(k)*(rb(k-1,iCell)*(qtot(k-1,iCell)) + &
+ rr(k-1,iCell)*(1.+qtot(k-1,iCell))) ))
+!SHP-old version
+! ( fzm(k)*rr(k,iCell) + fzm(k)*(rb(k,iCell)+rr(k,iCell))*qtot(k,iCell) &
+! +fzp(k)*rr(k-1,iCell) + fzp(k)*(rb(k-1,iCell)+rr(k-1,iCell))*qtot(k-1,iCell) ))
+
+
! - gravity*(fzm(k)*rb(k,iCell)+fzp(k)*rb(k-1,iCell)) ) &
! - gravity*( fzm(k)*(rr(k,iCell)+rb(k,iCell)) + fzp(k)*(rr(k-1,iCell)+rb(k-1,iCell)) )
-
-
! - cqw(k,iCell)*rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
! - gravity*( fzm(k)*rr(k,iCell)+fzp(k)*rr(k-1,iCell) &
! +(1.-cqw(k,iCell))*(fzm(k)*rb(k,iCell)+fzp(k)*rb(k-1,iCell)))
-
-
! WCS version - cqw(k,iCell)*rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
! - gravity*0.5*(rr(k,iCell)+rr(k-1,iCell)+(1.-cqw(k,iCell))*(rb(k,iCell)+rb(k-1,iCell)))
@@ -2587,7 +2630,7 @@
if ( v_mom_eddy_visc2 > 0.0 ) then
do iCell = 1, grid % nCellsSolve
- do k=2,nVertLevels-1
+ 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))*( &
(w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
-(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
@@ -2831,7 +2874,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) + rho(k,iCell)*h_diabatic(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + h_diabatic(k,iCell)
end do
end do
@@ -2873,8 +2916,8 @@
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))-(theta(k ,iCell)-t_init(k)))/(zp-z0) &
- -((theta(k ,iCell)-t_init(k))-(theta(k-1,iCell)-t_init(k-1)))/(z0-zm) )/(0.5*(zp-zm))
+ ((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))
end do
end do
@@ -3179,12 +3222,6 @@
enddo
enddo
- do i=1,grid%nCellsSolve
- do k=1,grid % nVertLevels + 1
- grid % rw % array (k,i) = 0.
- enddo
- enddo
-
end subroutine init_coupled_diagnostics
! ------------------------
@@ -3231,8 +3268,6 @@
do k = 1, grid % nVertLevels
- grid % rt_diabatic_tend % array(k,iCell) = state_new % theta % array(k,iCell)
-
state_new % theta % array(k,iCell) = t(k)*(1.+1.61*qv(k))
grid % rt_diabatic_tend % array(k,iCell) = state_new % rho % array(k,iCell) * &
(state_new % theta % array(k,iCell) - grid % rt_diabatic_tend % array(k,iCell))/dt
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-09-24 18:32:31 UTC (rev 510)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-09-24 22:57:36 UTC (rev 511)
@@ -7,6 +7,12 @@
type (domain_type), intent(inout) :: domain
+!
+! Note: initialize_advection_rk() and initialize_deformation_weights()
+! are called from within setup_nhyd_test_case(); after initialization
+! is migrated to a separate program/processing step, we can uncomment
+! calls to these routines in in mpas_init() if necessary
+!
call setup_nhyd_test_case(domain)
end subroutine mpas_setup_test_case
@@ -15,7 +21,7 @@
subroutine mpas_init(block, mesh, dt)
use grid_types
- use advection
+! use advection
use time_integration
implicit none
@@ -28,8 +34,12 @@
! call compute_state_diagnostics(block % state % time_levs(1) % state, mesh)
call init_coupled_diagnostics( block % state % time_levs(1) % state, mesh)
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh) ! ok for nonhydrostatic model
- call initialize_advection_rk(mesh)
- call initialize_deformation_weights(mesh)
+!
+! Note: The following initialization calls have been moved to mpas_setup_test_case()
+! since values computed by these routines are needed to produce initial fields
+!
+! call initialize_advection_rk(mesh)
+! call initialize_deformation_weights(mesh)
end subroutine mpas_init
</font>
</pre>