<p><b>laura@ucar.edu</b> 2010-10-07 12:50:32 -0600 (Thu, 07 Oct 2010)</p><p>reconciled with module_time_integration.F from /branches/atmos_nonhydrostatic version 520, except for stuff related to the reorganized registry, and stuff related to physics<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-10-07 18:48:01 UTC (rev 523)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-10-07 18:50:32 UTC (rev 524)
@@ -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 % time_levs(1) % state, config_test_case)
+ write(0,*) ' returned from test case setup '
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % 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
@@ -89,20 +108,25 @@
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 :: rh_max = 0.4 ! Maximum relative humidity
+ real (kind=RKIND), parameter :: rh_max = 0.65 ! 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
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 :: rh
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 +157,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,9 +180,19 @@
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
+ nCellsSolve = grid % nCellsSolve
zgrid => grid % zgrid % array
rdzw => grid % rdzw % array
@@ -183,9 +218,14 @@
rr => state % rho_p % array
t => state % theta % array
rt => grid % rtheta_p % array
+ rh => state % rh % array
scalars => state % scalars % array
+ call initialize_advection_rk(grid)
+ call initialize_deformation_weights(grid)
+
+ rh(:,:) = 0.
scalars(:,:,:) = 0.
xnutr = 0.
@@ -457,14 +497,14 @@
enddo
end if
+ 201 format(3i6,8(1x,e15.8))
+ 202 format(2i6,10(1x,e15.8))
#ifdef DO_PHYSICS
-!
!... initialization for moisture (following initialization in hydrostatic core):
! note that the initial profile of relative humidity is different from that
! in the hydrostatic core since the grid is different.
-
- 201 format(3i6,8(1x,e15.8))
- 202 format(2i6,8(1x,e15.8))
+! write(0,*)
+! write(0,*) 'RELATIVE HUMIDITY:'
do k = 1, nz1
ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
if(ztemp .gt. 5000.) then
@@ -475,7 +515,10 @@
rel_hum(k,i) = (1.-0.75*(ztemp/zt)**1.25)
endif
rel_hum(k,i) = min(rh_max,rel_hum(k,i))
+! if(ztemp .lt. 0.0) write(0,202) i,k,ztemp,rel_hum(k,i)
enddo
+#else
+ rel_hum(:,i) = 0.0
#endif
!
! iterations to converge temperature as a function of pressure
@@ -507,10 +550,8 @@
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.
#ifdef DO_PHYSICS
-!
!... computes the water vapor mixing ratio (following initialization in hydrostatic core):
! note that the initial profile of water vapor is different from that in the hydrostatic
! core since temperature and relative humidity are different.
@@ -520,11 +561,14 @@
es = 1000.*0.6112*exp(21.8745584*(tt(k)-273.16)/(tt(k)-7.66))
end if
qvs = (287.04/461.6)*es/(ptemp-es)
- scalars(index_qv,k,i) = rel_hum(k,i)*qvs
-! write(6,201) itr,i,k,ztemp,zgrid(k,i),zgrid(k+1,i),ptemp,tt(k),qvs,rel_hum(k,i),scalars(index_qv,k,i)
+ rh(k,i) = rel_hum(k,i)
+ qv(k,i) = amin1(0.014,rel_hum(k,i)*qvs)
+ scalars(index_qv,k,i) = amin1(0.014,rel_hum(k,i)*qvs)
+#else
+ qv(k,i) = 0.
+ scalars(index_qv,k,i) = 0.
#endif
end do
-
! do k=2,nz1
! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,i)))
! end do
@@ -557,6 +601,7 @@
do k=1,nz1
p (k,i) = ((ppb(k,i)+pp(k,i))/p0)**(rgas/cp)
t (k,i) = tt(k)/p(k,i)
+! t (k,i) = tt(k)*(1.+1.61*scalars(index_qv,k,i))/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)
end do
@@ -615,6 +660,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
!
@@ -635,12 +688,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)
!
@@ -733,50 +863,37 @@
type (grid_state), 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
-
!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
- 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
!
! Scale all distances
@@ -839,12 +956,12 @@
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)
@@ -866,7 +983,7 @@
zt = 20000.
dz = zt/float(nz1)
- write(0,*) ' dz = ',dz
+! write(0,*) ' dz = ',dz
write(0,*) ' hx computation complete '
do k=1,nz
@@ -893,7 +1010,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)
@@ -956,15 +1073,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
!
@@ -973,22 +1081,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))
@@ -1001,6 +1105,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
@@ -1029,72 +1137,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)) &
@@ -1102,67 +1265,574 @@
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 (grid_meta), intent(inout) :: grid
+ type (grid_state), 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
+
+ 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
+
+ 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
@@ -1176,33 +1846,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
@@ -1226,16 +1955,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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
</font>
</pre>