<p><b>duda</b> 2012-07-25 18:08:47 -0600 (Wed, 25 Jul 2012)</p><p>BRANCH COMMIT<br>
<br>
Update initialization and model to use the scaled planet radius that is stored <br>
in grid % sphere_radius, which is set when initial conditions are generated.<br>
<br>
<br>
M src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
M src/core_nhyd_atmos/mpas_atm_advection.F<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-07-25 03:41:33 UTC (rev 2055)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-07-26 00:08:47 UTC (rev 2056)
@@ -263,22 +263,22 @@
logical, parameter :: moisture = .true.
!
- ! Scale all distances and areas from a unit sphere to one with radius a
+ ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
!
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
+ grid % xCell % array = grid % xCell % array * grid % sphere_radius
+ grid % yCell % array = grid % yCell % array * grid % sphere_radius
+ grid % zCell % array = grid % zCell % array * grid % sphere_radius
+ grid % xVertex % array = grid % xVertex % array * grid % sphere_radius
+ grid % yVertex % array = grid % yVertex % array * grid % sphere_radius
+ grid % zVertex % array = grid % zVertex % array * grid % sphere_radius
+ grid % xEdge % array = grid % xEdge % array * grid % sphere_radius
+ grid % yEdge % array = grid % yEdge % array * grid % sphere_radius
+ grid % zEdge % array = grid % zEdge % array * grid % sphere_radius
+ grid % dvEdge % array = grid % dvEdge % array * grid % sphere_radius
+ grid % dcEdge % array = grid % dcEdge % array * grid % sphere_radius
+ grid % areaCell % array = grid % areaCell % array * grid % sphere_radius**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * grid % sphere_radius**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * grid % sphere_radius**2.0
weightsOnEdge => grid % weightsOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
@@ -344,7 +344,7 @@
znut = eta_t
etavs = (1.-0.252)*pii/2.
- r_earth = a
+ r_earth = grid % sphere_radius
omega_e = omega
p0 = 1.e+05
@@ -589,7 +589,7 @@
end do
call init_atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
+ cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat,grid%sphere_radius)
end if
@@ -746,24 +746,24 @@
lat2 = grid%latVertex%array(vtx2)
iCell1 = grid % cellsOnEdge % array(1,iEdge)
iCell2 = grid % cellsOnEdge % array(2,iEdge)
- flux = (0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
+ flux = (0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1))) * grid % sphere_radius / grid % dvEdge % array(iEdge)
if (config_test_case == 2) then
r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &
lat_pert, lon_pert, 1.0_RKIND)/(pert_radius)
- u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
+ u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1) * grid % sphere_radius / grid % dvEdge % array(iEdge)
else if (config_test_case == 3) then
lon_Edge = grid % lonEdge % array(iEdge)
u_pert = u_perturbation*cos(k_x*(lon_Edge - lon_pert)) &
- *(0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
+ *(0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1))) * grid % sphere_radius / grid % dvEdge % array(iEdge)
else
u_pert = 0.0
end if
if (rebalance) then
- call init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),a,u0,nz1,nlat)
+ call init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),grid%sphere_radius,u0,nz1,nlat)
do k=1,grid % nVertLevels
fluxk = u0*flux_zonal(k)/(0.5*(rb(k,iCell1)+rb(k,iCell2)+rr(k,iCell1)+rr(k,iCell2)))
state % u % array(k,iEdge) = fluxk + u_pert
@@ -923,6 +923,7 @@
end subroutine init_atm_test_case_jw
+
subroutine init_atm_calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
implicit none
@@ -972,9 +973,11 @@
end subroutine init_atm_calc_flux_zonal
+
+
!SHP-balance
subroutine init_atm_recompute_geostrophic_wind(u_2d,rho_2d,pp_2d,qv_2d,lat_2d,zz_2d,zx_2d, &
- cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat)
+ cf1,cf2,cf3,fzm,fzp,rdzw,nz1,nlat,dlat,rad)
implicit none
integer, intent(in) :: nz1,nlat
@@ -983,7 +986,7 @@
real (kind=RKIND), dimension(nz1,nlat-1), intent(in) :: zx_2d
real (kind=RKIND), dimension(nlat), intent(in) :: lat_2d
real (kind=RKIND), dimension(nz1), intent(in) :: fzm, fzp, rdzw
- real (kind=RKIND), intent(in) :: cf1, cf2, cf3, dlat
+ real (kind=RKIND), intent(in) :: cf1, cf2, cf3, dlat, rad
!local variable
real (kind=RKIND), dimension(nz1,nlat-1) :: pgrad, ru, u
@@ -996,7 +999,7 @@
real (kind=RKIND) :: rdx, qtot, r_earth, phi
integer :: k,i, itr
- r_earth = a
+ r_earth = rad
omega_e = omega
rdx = 1./(dlat*r_earth)
@@ -1067,10 +1070,9 @@
u_2d(k,i) = (3.*u_2d(k,i-1)-u_2d(k,i-2))*.5
end do
-
end subroutine init_atm_recompute_geostrophic_wind
-!----------------------------------------------------------------------------------------------------------
+
subroutine init_atm_test_case_squall_line(dminfo, grid, state, diag, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup squall line and supercell test case
@@ -2395,7 +2397,7 @@
etavs = (1.-0.252)*pii/2.
rcv = rgas/(cp-rgas)
- r_earth = a
+ r_earth = grid % sphere_radius
omega_e = omega
p0 = 1.e+05
@@ -2405,25 +2407,25 @@
!
- ! Scale all distances and areas from a unit sphere to one with radius a
+ ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
!
if (config_static_interp) then
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
+ grid % xCell % array = grid % xCell % array * r_earth
+ grid % yCell % array = grid % yCell % array * r_earth
+ grid % zCell % array = grid % zCell % array * r_earth
+ grid % xVertex % array = grid % xVertex % array * r_earth
+ grid % yVertex % array = grid % yVertex % array * r_earth
+ grid % zVertex % array = grid % zVertex % array * r_earth
+ grid % xEdge % array = grid % xEdge % array * r_earth
+ grid % yEdge % array = grid % yEdge % array * r_earth
+ grid % zEdge % array = grid % zEdge % array * r_earth
+ grid % dvEdge % array = grid % dvEdge % array * r_earth
+ grid % dcEdge % array = grid % dcEdge % array * r_earth
+ grid % areaCell % array = grid % areaCell % array * r_earth**2.0
+ grid % areaTriangle % array = grid % areaTriangle % array * r_earth**2.0
+ grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * r_earth**2.0
scalars(:,:,:) = 0.
@@ -4362,6 +4364,7 @@
end subroutine init_atm_test_case_gfs
+
subroutine init_atm_test_case_sfc(domain, dminfo, grid, fg, state, diag, test_case, parinfo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Real-data test case using SST data
@@ -4656,10 +4659,8 @@
!
! Scale all distances
!
+ a_scale = grid % sphere_radius
-! a_scale = 1.0
- a_scale = a
-
grid % xCell % array = grid % xCell % array * a_scale
grid % yCell % array = grid % yCell % array * a_scale
grid % zCell % array = grid % zCell % array * a_scale
@@ -4804,7 +4805,7 @@
grid % cf2 % scalar = cf2
grid % cf3 % scalar = cf3
- write(0,*) 'EARTH RADIUS = ', a
+ write(0,*) 'EARTH RADIUS = ', grid % sphere_radius
! setting for terrain
@@ -4816,18 +4817,18 @@
xi = grid % lonCell % array(iCell)
yi = grid % latCell % array(iCell)
- xc = sphere_distance(yi, xi, yi, 0., a)
- yc = sphere_distance(yi, xi, 0., xi, a)
- xac = sphere_distance(yi, xa /a, yi, 0., a)
- xlac = sphere_distance(yi, xla/a, yi, 0., a)
+ xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
+ yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
+ xac = sphere_distance(yi, xa /grid % sphere_radius, yi, 0., grid % sphere_radius)
+ xlac = sphere_distance(yi, xla/grid % sphere_radius, yi, 0., grid % sphere_radius)
- ri = sphere_distance(yi, xi, 0., 0., a)
+ ri = sphere_distance(yi, xi, 0., 0., grid % sphere_radius)
! Circular mountain with Schar mtn cross section
hx(1,iCell) = hm*exp(-(ri/xa)**2)*cos(pii*ri/xla)**2
! Ridge mountain with Schar mtn cross section
-! hx(1,iCell) = hm*exp(-(xc/xac)**2)*cos(pii*xc/xlac)**2*cos(yc/a)
+! hx(1,iCell) = hm*exp(-(xc/xac)**2)*cos(pii*xc/xlac)**2*cos(yc/grid % sphere_radius)
hx(nz,iCell) = zt
@@ -4968,7 +4969,7 @@
!
allocate(psiVertex(grid % nVertices))
do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * um * ( &
+ psiVertex(iVtx) = -grid % sphere_radius * um * ( &
sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
)
@@ -5325,10 +5326,8 @@
!
! Scale all distances
!
+ a_scale = grid % sphere_radius
-! a_scale = 1.0
- a_scale = a
-
grid % xCell % array = grid % xCell % array * a_scale
grid % yCell % array = grid % yCell % array * a_scale
grid % zCell % array = grid % zCell % array * a_scale
@@ -5468,7 +5467,7 @@
write(6,*) k,zw(k), ah(k)
end do
- write(0,*) 'EARTH RADIUS = ', a
+ write(0,*) 'EARTH RADIUS = ', grid % sphere_radius
! for hx computation
r1m = .75*pii
@@ -5478,7 +5477,7 @@
xa = pii/16. ! for specifying mtn with in degrees
- xa = pii*a/16. ! corresponds to ~11 grid intervals across entire mtn with 2 deg res
+ xa = pii*grid%sphere_radius/16. ! corresponds to ~11 grid intervals across entire mtn with 2 deg res
do iCell=1,grid % nCells
@@ -5499,16 +5498,16 @@
xi = grid % lonCell % array(iCell)
yi = grid % latCell % array(iCell)
- xc = sphere_distance(yi, xi, yi, 0., a)
- yc = sphere_distance(yi, xi, 0., xi, a)
+ xc = sphere_distance(yi, xi, yi, 0., grid % sphere_radius)
+ yc = sphere_distance(yi, xi, 0., xi, grid % sphere_radius)
if(abs(xc).ge.xa) then
! if(abs(xi).ge.xa.and.abs(2.*pii-xi).ge.xa) then
hx(1,iCell) = 0.
else
! for mtn ridge with uniform width in km
- hx(1,iCell) = hm*cos(.5*pii*xc/xa)**2*cos(yc/a)
+ hx(1,iCell) = hm*cos(.5*pii*xc/xa)**2*cos(yc/grid % sphere_radius)
! for mtn ridge with uniform width in degrees
-! hx(1,iCell) = hm*cos(.5*pii*xi/xa)**2*cos(yc/a)
+! hx(1,iCell) = hm*cos(.5*pii*xi/xa)**2*cos(yc/grid % sphere_radius)
end if
hx(:,iCell) = hx(1,iCell)
Modified: branches/dcmip/src/core_nhyd_atmos/mpas_atm_advection.F
===================================================================
--- branches/dcmip/src/core_nhyd_atmos/mpas_atm_advection.F 2012-07-25 03:41:33 UTC (rev 2055)
+++ branches/dcmip/src/core_nhyd_atmos/mpas_atm_advection.F 2012-07-26 00:08:47 UTC (rev 2056)
@@ -111,9 +111,9 @@
do i=1,n
advCells(i+1,iCell) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1,iCell))/a
- yc(i) = grid % yCell % array(advCells(i+1,iCell))/a
- zc(i) = grid % zCell % array(advCells(i+1,iCell))/a
+ xc(i) = grid % xCell % array(advCells(i+1,iCell))/grid%sphere_radius
+ yc(i) = grid % yCell % array(advCells(i+1,iCell))/grid%sphere_radius
+ zc(i) = grid % zCell % array(advCells(i+1,iCell))/grid%sphere_radius
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
@@ -131,8 +131,8 @@
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
+ dl_sphere(i) = grid%sphere_radius*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
end do
length_scale = 1.
@@ -262,12 +262,12 @@
if (ip1 > n-1) ip1 = 1
iEdge = grid % EdgesOnCell % array (i,iCell)
- xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
+ xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/grid%sphere_radius
+ xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
+ yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
+ zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/grid%sphere_radius
if ( grid % on_a_sphere ) then
call arc_bisect( xv1, yv1, zv1, &
@@ -825,16 +825,16 @@
! compute poynomial fit for this cell if all needed neighbors exist
if (grid % on_a_sphere) then
- xc(1) = grid % xCell % array(iCell)/a
- yc(1) = grid % yCell % array(iCell)/a
- zc(1) = grid % zCell % array(iCell)/a
+ xc(1) = grid % xCell % array(iCell)/grid%sphere_radius
+ yc(1) = grid % yCell % array(iCell)/grid%sphere_radius
+ zc(1) = grid % zCell % array(iCell)/grid%sphere_radius
do i=2,n
iv = grid % verticesOnCell % array(i-1,iCell)
- xc(i) = grid % xVertex % array(iv)/a
- yc(i) = grid % yVertex % array(iv)/a
- zc(i) = grid % zVertex % array(iv)/a
+ xc(i) = grid % xVertex % array(iv)/grid%sphere_radius
+ yc(i) = grid % yVertex % array(iv)/grid%sphere_radius
+ zc(i) = grid % zVertex % array(iv)/grid%sphere_radius
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
@@ -852,8 +852,8 @@
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
+ dl_sphere(i) = grid%sphere_radius*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
end do
length_scale = 1.
Modified: branches/dcmip/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/dcmip/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-07-25 03:41:33 UTC (rev 2055)
+++ branches/dcmip/src/core_nhyd_atmos/mpas_atm_time_integration.F 2012-07-26 00:08:47 UTC (rev 2056)
@@ -1861,7 +1861,7 @@
!SHP-curvature
logical, parameter :: curvature = .true.
!real (kind=RKIND), parameter :: omega_e = 7.29212e-05*100.
- !real (kind=RKIND) :: r_earth
+ real (kind=RKIND) :: r_earth
real (kind=RKIND), dimension(:,:), pointer :: ur_cell, vr_cell
real (kind=RKIND), parameter :: c_s = 0.125
@@ -1886,7 +1886,7 @@
!-----------
!SHP-curvature
- !r_earth = a
+ r_earth = grid % sphere_radius
ur_cell => diag % uReconstructZonal % array
vr_cell => diag % uReconstructMeridional % array
@@ -2153,7 +2153,7 @@
- 2.*omega*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
*rho_edge(k,iEdge)*.25*(w(k,cell1)+w(k+1,cell1)+w(k,cell2)+w(k+1,cell2)) &
- u(k,iEdge)*.25*(w(k+1,cell1)+w(k,cell1)+w(k,cell2)+w(k+1,cell2)) &
- *rho_edge(k,iEdge)/a
+ *rho_edge(k,iEdge)/r_earth
!old-err.
!tend_u(k,iEdge) = tend_u(k,iEdge) &
! - 2.*omega_e*cos(grid % angleEdge % array(iEdge))*cos(grid % latEdge % array(iEdge)) &
@@ -2512,7 +2512,7 @@
do k=2,nVertLevels
tend_w(k,iCell) = tend_w(k,iCell) + (rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k))* &
( (fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell))**2. &
- +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/a &
+ +(fzm(k)*vr_cell(k,iCell)+fzp(k)*vr_cell(k-1,iCell))**2. )/r_earth &
+ 2.*omega*cos(grid % latCell % array(iCell)) &
*(fzm(k)*ur_cell(k,iCell)+fzp(k)*ur_cell(k-1,iCell)) &
*(rho_zz(k,iCell)*fzm(k)+rho_zz(k-1,iCell)*fzp(k))
</font>
</pre>