<p><b>duda</b> 2013-04-17 18:25:31 -0600 (Wed, 17 Apr 2013)</p><p>BRANCH COMMIT<br>
<br>
Cleanup of grid dereferences in mpas_atm_time_integration.F. Primarily, change<br>
lines like<br>
<br>
do iCell=1,grid%nCells<br>
<br>
to<br>
<br>
do iCell=1,nCells<br>
<br>
No impact on results.<br>
<br>
<br>
M src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F 2013-04-17 21:03:20 UTC (rev 2765)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F 2013-04-18 00:25:31 UTC (rev 2766)
@@ -354,8 +354,8 @@
do k = 1, block % mesh % nVertLevels
scalar_min = min(scalar_min, block % state % time_levs(2) % state % w % array(k,iCell))
scalar_max = max(scalar_max, block % state % time_levs(2) % state % w % array(k,iCell))
- enddo
- enddo
+ end do
+ end do
call mpas_dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
call mpas_dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
write(0,*) 'global min, max w ',global_scalar_min, global_scalar_max
@@ -366,8 +366,8 @@
do k = 1, block % mesh % nVertLevels
scalar_min = min(scalar_min, block % state % time_levs(2) % state % u % array(k,iEdge))
scalar_max = max(scalar_max, block % state % time_levs(2) % state % u % array(k,iEdge))
- enddo
- enddo
+ end do
+ end do
call mpas_dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
call mpas_dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
write(0,*) 'global min, max u ',global_scalar_min, global_scalar_max
@@ -379,8 +379,8 @@
do k = 1, block % mesh % nVertLevels
scalar_min = min(scalar_min, block % state % time_levs(2) % state % scalars % array(iScalar,k,iCell))
scalar_max = max(scalar_max, block % state % time_levs(2) % state % scalars % array(iScalar,k,iCell))
- enddo
- enddo
+ end do
+ end do
call mpas_dmpar_min_real(domain%dminfo, scalar_min, global_scalar_min)
call mpas_dmpar_max_real(domain%dminfo, scalar_max, global_scalar_max)
write(0,102) iScalar,global_scalar_min,global_scalar_max
@@ -429,12 +429,15 @@
integer :: iEdge, iCell, k, cell1, cell2, iq
integer :: nCells, nEdges, nVertLevels, nCellsSolve
real (kind=RKIND) :: qtot
+ integer, dimension(:,:), pointer :: cellsOnEdge
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
nCellsSolve = grid % nCellsSolve
+ cellsOnEdge => grid % cellsOnEdge % array
+
do iCell = 1, nCellsSolve
do k = 2, nVertLevels
qtot = 0.
@@ -446,8 +449,8 @@
end do
do iEdge = 1, nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k = 1, nVertLevels
qtot = 0.
@@ -481,7 +484,7 @@
real (kind=RKIND), intent(in) :: dts
- integer :: i, k, iq
+ integer :: iCell, k, iq
integer :: nCells, nVertLevels, nCellsSolve
real (kind=RKIND), dimension(:,:), pointer :: zz, cqw, p, t, rb, rtb, pb, rt
@@ -530,52 +533,52 @@
cofrz(k) = dtseps*rdzw(k)
end do
- do i = 1, nCellsSolve ! we only need to do cells we are solving for, not halo cells
+ do iCell = 1, nCellsSolve ! we only need to do cells we are solving for, not halo cells
do k=2,nVertLevels
- cofwr(k,i) =.5*dtseps*gravity*(fzm(k)*zz(k,i)+fzp(k)*zz(k-1,i))
+ cofwr(k,iCell) =.5*dtseps*gravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))
end do
- coftz(1,i) = 0.0
+ coftz(1,iCell) = 0.0
do k=2,nVertLevels
- cofwz(k,i) = dtseps*c2*(fzm(k)*zz(k,i)+fzp(k)*zz(k-1,i)) &
- *rdzu(k)*cqw(k,i)*(fzm(k)*p (k,i)+fzp(k)*p (k-1,i))
- coftz(k,i) = dtseps* (fzm(k)*t (k,i)+fzp(k)*t (k-1,i))
+ cofwz(k,iCell) = dtseps*c2*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell)) &
+ *rdzu(k)*cqw(k,iCell)*(fzm(k)*p (k,iCell)+fzp(k)*p (k-1,iCell))
+ coftz(k,iCell) = dtseps* (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell))
end do
- coftz(nVertLevels+1,i) = 0.0
+ coftz(nVertLevels+1,iCell) = 0.0
do k=1,nVertLevels
qtot = 0.
do iq = s % moist_start, s % moist_end
- qtot = qtot + s % scalars % array (iq, k, i)
+ qtot = qtot + s % scalars % array (iq, k, iCell)
end do
- cofwt(k,i) = .5*dtseps*rcv*zz(k,i)*gravity*rb(k,i)/(1.+qtot) &
- *p(k,i)/((rtb(k,i)+rt(k,i))*pb(k,i))
+ cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtot) &
+ *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell))
end do
- a_tri(1,i) = 0. ! note, this value is never used
+ a_tri(1,iCell) = 0. ! note, this value is never used
b_tri(1) = 1. ! note, this value is never used
c_tri(1) = 0. ! note, this value is never used
- gamma_tri(1,i) = 0.
- alpha_tri(1,i) = 0. ! note, this value is never used
+ gamma_tri(1,iCell) = 0.
+ alpha_tri(1,iCell) = 0. ! note, this value is never used
do k=2,nVertLevels
- a_tri(k,i) = -cofwz(k ,i)* coftz(k-1,i)*rdzw(k-1)*zz(k-1,i) &
- +cofwr(k ,i)* cofrz(k-1 ) &
- -cofwt(k-1,i)* coftz(k-1,i)*rdzw(k-1)
+ a_tri(k,iCell) = -cofwz(k ,iCell)* coftz(k-1,iCell)*rdzw(k-1)*zz(k-1,iCell) &
+ +cofwr(k ,iCell)* cofrz(k-1 ) &
+ -cofwt(k-1,iCell)* coftz(k-1,iCell)*rdzw(k-1)
b_tri(k) = 1. &
- +cofwz(k ,i)*(coftz(k ,i)*rdzw(k )*zz(k ,i) &
- +coftz(k ,i)*rdzw(k-1)*zz(k-1,i)) &
- -coftz(k ,i)*(cofwt(k ,i)*rdzw(k ) &
- -cofwt(k-1,i)*rdzw(k-1)) &
- +cofwr(k, i)*(cofrz(k )-cofrz(k-1))
- c_tri(k) = -cofwz(k ,i)* coftz(k+1,i)*rdzw(k )*zz(k ,i) &
- -cofwr(k ,i)* cofrz(k ) &
- +cofwt(k ,i)* coftz(k+1,i)*rdzw(k )
+ +cofwz(k ,iCell)*(coftz(k ,iCell)*rdzw(k )*zz(k ,iCell) &
+ +coftz(k ,iCell)*rdzw(k-1)*zz(k-1,iCell)) &
+ -coftz(k ,iCell)*(cofwt(k ,iCell)*rdzw(k ) &
+ -cofwt(k-1,iCell)*rdzw(k-1)) &
+ +cofwr(k, iCell)*(cofrz(k )-cofrz(k-1))
+ c_tri(k) = -cofwz(k ,iCell)* coftz(k+1,iCell)*rdzw(k )*zz(k ,iCell) &
+ -cofwr(k ,iCell)* cofrz(k ) &
+ +cofwt(k ,iCell)* coftz(k+1,iCell)*rdzw(k )
end do
do k=2,nVertLevels
- alpha_tri(k,i) = 1./(b_tri(k)-a_tri(k,i)*gamma_tri(k-1,i))
- gamma_tri(k,i) = c_tri(k)*alpha_tri(k,i)
+ alpha_tri(k,iCell) = 1./(b_tri(k)-a_tri(k,iCell)*gamma_tri(k-1,iCell))
+ gamma_tri(k,iCell) = c_tri(k)*alpha_tri(k,iCell)
end do
end do ! loop over cells
@@ -599,6 +602,7 @@
type (mesh_type) :: grid
integer :: iCell, iEdge, k, cell1, cell2, coef_3rd_order
+ integer :: nCellsSolve, nCells, nVertLevels, nEdges
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
real (kind=RKIND) :: flux
@@ -606,6 +610,11 @@
coef_3rd_order = config_coef_3rd_order
if (config_theta_adv_order /=3) coef_3rd_order = 0
+ nCellsSolve = grid % nCellsSolve
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
fzm => grid % fzm % array
fzp => grid % fzp % array
dvEdge => grid % dvEdge % array
@@ -624,8 +633,8 @@
! we solve for omega instead of w (see Klemp et al MWR 2007),
! so here we change the w_p tendency to an omega_p tendency
- do iCell = 1, grid % nCellsSolve
- do k = 2, grid % nVertLevels
+ do iCell = 1, nCellsSolve
+ do k = 2, nVertLevels
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)
@@ -635,12 +644,12 @@
! here we need to compute the omega tendency in a manner consistent with our diagnosis of omega.
! this requires us to use the same flux divergence as is used in the theta eqn - see Klemp et al MWR 2003.
- do iEdge = 1,grid % nEdges
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k = 2, grid%nVertLevels
+ do k = 2, 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) &
+ (grid % zb % array(k,2,iEdge) + coef_3rd_order*sign(1.0_RKIND,tend % u % array(k,iEdge))*grid %zb3 % array(k,2,iEdge))*flux &
@@ -687,6 +696,7 @@
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, rdzw, dcEdge, AreaCell, cofrz, dvEdge
real (kind=RKIND), dimension(:,:), pointer :: cpr, cpl, pzp, pzm
+ integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND) :: smdiv, c2, rcv
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
@@ -706,6 +716,7 @@
logical :: newpx
!--
+ cellsOnEdge => grid % cellsOnEdge % array
rho_zz => s % rho_zz % array
theta_m => s % theta_m % array
@@ -789,11 +800,11 @@
do iEdge = 1, nEdges
- cell1 = grid % cellsOnEdge % array (1,iEdge)
- cell2 = grid % cellsOnEdge % array (2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
! update edges for block-owned cells
- if (cell1 <= grid % nCellsSolve .or. cell2 <= grid % nCellsSolve ) then
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
if (newpx) then
@@ -839,7 +850,7 @@
+pzp(k,cell1)*(zz(k+2,cell1)*rtheta_pp_old(k+2,cell1) &
-zz(k ,cell1)*rtheta_pp_old(k ,cell1)))
- do k=2,grid % nVertLevels-1
+ do k=2,nVertLevels-1
dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
*(pzp(k,cell2)*(zz(k+1,cell2)*rtheta_pp_old(k+1,cell2) &
-zz(k ,cell2)*rtheta_pp_old(k ,cell2)) &
@@ -851,7 +862,7 @@
-zz(k-1,cell1)*rtheta_pp_old(k-1,cell1)))
end do
- k = grid % nVertLevels
+ k = nVertLevels
dpzx(k) = .25*(zx(k,iEdge)+zx(k+1,iEdge)) &
*(pzm(k,cell2)*(zz(k ,cell2)*rtheta_pp_old(k ,cell2) &
-zz(k-1,cell2)*rtheta_pp_old(k-1,cell2)) &
@@ -986,7 +997,7 @@
rho_pp, rho_zz, rho_base, ruAvg, ru_save, ru_p, u, ru, &
exner, exner_base, rtheta_base, pressure_p, &
zz, theta_m, pressure_b, qvapor
- real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell
+ real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -1035,8 +1046,8 @@
fzm => grid % fzm % array
fzp => grid % fzp % array
dvEdge => grid % dvEdge % array
- AreaCell => grid % AreaCell % array
- CellsOnEdge => grid % CellsOnEdge % array
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
nVertLevels = grid % nVertLevels
nCells = grid % nCells
@@ -1105,8 +1116,8 @@
do iEdge = 1, nEdges
- cell1 = CellsOnEdge(1,iEdge)
- cell2 = CellsOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
do k = 1, nVertLevels
ruAvg(k,iEdge) = ru(k,iEdge) + (ruAvg(k,iEdge) / float(ns))
@@ -1182,7 +1193,7 @@
real (kind=RKIND), dimension( s_old % num_scalars, grid % nVertLevels ) :: flux_arr
real (kind=RKIND), dimension( s_old % num_scalars, grid % nVertLevels + 1 ) :: wdtn
- integer :: nVertLevels
+ integer :: nCellsSolve, nEdges, nVertLevels
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
real (kind=RKIND) :: coef_3rd_order
@@ -1228,6 +1239,8 @@
adv_coefs => grid % adv_coefs % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
h_theta_eddy_visc2 = config_h_theta_eddy_visc2
@@ -1246,10 +1259,10 @@
!
! horizontal flux divergence, accumulate in scalar_tend
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then ! only for owned cells
! flux_arr stores the value of the scalar at the edge.
! a better name perhaps would be scalarEdge
@@ -1257,7 +1270,7 @@
flux_arr(:,:) = 0.
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge)
do iScalar=1,s_old % num_scalars
flux_arr(iScalar,k) = flux_arr(iScalar,k) + scalar_weight* scalar_new(iScalar,k,iCell)
@@ -1268,7 +1281,7 @@
! here we add the horizontal flux divergence into the scalar tendency.
! note that the scalar tendency is modified.
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
do iScalar=1,s_old % num_scalars
scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) &
- uhAvg(k,iEdge)*flux_arr(iScalar,k)/areaCell(cell1)
@@ -1287,17 +1300,17 @@
! zero fluxes at top and bottom
wdtn(:,1) = 0.
- wdtn(:,grid % nVertLevels+1) = 0.
+ wdtn(:,nVertLevels+1) = 0.
if (config_scalar_vadv_order == 2) then
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
do k = 2, nVertLevels
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
end do
end do
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1307,12 +1320,12 @@
else if ( config_scalar_vadv_order == 3 ) then
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
k = 2
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- enddo
+ end do
do k=3,nVertLevels-1
do iScalar=1,s_old % num_scalars
@@ -1324,9 +1337,9 @@
k = nVertLevels
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- enddo
+ end do
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1337,12 +1350,12 @@
else if ( config_scalar_vadv_order == 4 ) then
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
k = 2
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- enddo
+ end do
do k=3,nVertLevels-1
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = flux4( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
@@ -1352,9 +1365,9 @@
k = nVertLevels
do iScalar=1,s_old % num_scalars
wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
- enddo
+ end do
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
do iScalar=1,s_old % num_scalars
scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
@@ -1430,9 +1443,10 @@
real (kind=RKIND), dimension( grid % nVertLevels, grid % nEdges ) :: flux_arr
real (kind=RKIND), dimension( grid % nVertLevels + 1, grid % nCells ) :: wdtn
- integer :: nVertLevels, num_scalars, icellmax, kmax
+ integer :: nCells, nCellsSolve, nEdges, nVertLevels, num_scalars, icellmax, kmax
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
+ integer, dimension(:), pointer :: nEdgesOnCell
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2
@@ -1474,11 +1488,15 @@
meshScalingDel2 => grid % meshScalingDel2 % array
meshScalingDel4 => grid % meshScalingDel4 % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
h_theta_eddy_visc2 = config_h_theta_eddy_visc2
@@ -1497,8 +1515,8 @@
! Note, however, that we enforce positive-definiteness in this update.
! The transport will maintain this positive definite solution and optionally, shape preservation (monotonicity).
- do iCell = 1,grid%nCellsSolve
- do k = 1, grid%nVertLevels
+ do iCell = 1, nCellsSolve
+ do k = 1, nVertLevels
do iScalar = 1,s_old%num_scalars
s_old % scalars % array(iScalar,k,iCell) = s_old % scalars % array(iScalar,k,iCell)+dt*scalar_tend(iScalar,k,iCell) / h_old(k,iCell)
scalar_tend(iScalar,k,iCell) = 0.
@@ -1521,8 +1539,8 @@
do iScalar = 1, s_old % num_scalars
write(0,*) ' mono transport for scalar ',iScalar
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
+ do iCell = 1, nCells
+ do k = 1, nVertLevels
scalar_old(k,iCell) = s_old % scalars % array(iScalar,k,iCell)
scalar_new(k,iCell) = s_new % scalars % array(iScalar,k,iCell)
end do
@@ -1531,22 +1549,22 @@
if (debug_print) then
scmin = scalar_old(1,1)
scmax = scalar_old(1,1)
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
+ do iCell = 1, nCells
+ do k=1, nVertLevels
scmin = min(scmin,scalar_old(k,iCell))
scmax = max(scmax,scalar_old(k,iCell))
- enddo
- enddo
+ end do
+ end do
write(0,*) ' scmin, scmin old in ',scmin,scmax
scmin = scalar_new(1,1)
scmax = scalar_new(1,1)
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
+ do iCell = 1, nCells
+ do k=1, nVertLevels
scmin = min(scmin,scalar_new(k,iCell))
scmax = max(scmax,scalar_new(k,iCell))
- enddo
- enddo
+ end do
+ end do
write(0,*) ' scmin, scmin new in ',scmin,scmax
end if
@@ -1556,7 +1574,7 @@
!
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
! zero flux at top and bottom
wdtn(1,iCell) = 0.
@@ -1586,7 +1604,7 @@
! pull s_min and s_max from the (horizontal) surrounding cells
- do i=1, grid % nEdgesOnCell % array(iCell)
+ do i=1, nEdgesOnCell(iCell)
do k=1, grid % nVertLevels
s_max(k,iCell) = max(s_max(k,iCell),scalar_old(k, grid % CellsOnCell % array(i,iCell)))
s_min(k,iCell) = min(s_min(k,iCell),scalar_old(k, grid % CellsOnCell % array(i,iCell)))
@@ -1599,16 +1617,16 @@
! horizontal flux divergence
flux_arr(:,:) = 0.
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then ! only for owned cells
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
scalar_weight = uhAvg(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge))
flux_arr(k,iEdge) = flux_arr(k,iEdge) + scalar_weight* scalar_new(k,iCell)
end do
@@ -1620,7 +1638,7 @@
! vertical flux divergence for upwind update, we will put upwind update into scalar_new, and put factor of dt in fluxes
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
k = 1
scalar_new(k,iCell) = scalar_old(k,iCell)*h_old(k,iCell)
@@ -1647,11 +1665,11 @@
! upwind flux computation
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
- do k=1,grid % nVertLevels
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then ! only for owned cells
+ do k=1, nVertLevels
flux_upwind = grid % dvEdge % array(iEdge) * dt * &
(max(0.0_RKIND,uhAvg(k,iEdge))*scalar_old(k,cell1) + min(0.0_RKIND,uhAvg(k,iEdge))*scalar_old(k,cell2))
flux_arr(k,iEdge) = dt*flux_arr(k,iEdge) - flux_upwind
@@ -1669,7 +1687,7 @@
! next, the limiter
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k = 1, nVertLevels
s_min_update = (scalar_new(k,iCell)+scale_arr(SCALE_OUT,k,iCell))/h_new(k,iCell)
s_max_update = (scalar_new(k,iCell)+scale_arr(SCALE_IN,k,iCell))/h_new(k,iCell)
@@ -1705,10 +1723,10 @@
!
! rescale the fluxes
!
- do iEdge = 1, grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then
+ do iEdge = 1, nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
do k = 1, nVertLevels
flux = flux_arr(k,iEdge)
flux = max(0.0_RKIND,flux) * min(scale_arr(SCALE_OUT,k,cell1), scale_arr(SCALE_IN, k,cell2)) &
@@ -1720,7 +1738,7 @@
! rescale the vertical flux
- do iCell=1,grid % nCells
+ do iCell=1,nCells
do k = 2, nVertLevels
flux = wdtn(k,iCell)
flux = max(0.0_RKIND,flux) * min(scale_arr(SCALE_OUT,k-1,iCell), scale_arr(SCALE_IN,k ,iCell)) &
@@ -1731,19 +1749,19 @@
!
! do the scalar update now that we have the fluxes
!
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
- do k=1,grid % nVertLevels
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then ! only for owned cells
+ do k=1,nVertLevels
scalar_new(k,cell1) = scalar_new(k,cell1) - flux_arr(k,iEdge)/areaCell(cell1)
scalar_new(k,cell2) = scalar_new(k,cell2) + flux_arr(k,iEdge)/areaCell(cell2)
end do
end if
end do
- do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevels
+ do iCell=1,nCellsSolve
+ do k=1,nVertLevels
scalar_new(k,iCell) = ( scalar_new(k,iCell) &
+ (-rdnw(k)*(wdtn(k+1,iCell)-wdtn(k,iCell)) ) )/h_new(k,iCell)
end do
@@ -1753,8 +1771,8 @@
scmin = scalar_new(1,1)
scmax = scalar_new(1,1)
- do iCell = 1, grid%nCellsSolve
- do k=1, grid%nVertLevels
+ do iCell = 1, nCellsSolve
+ do k=1, nVertLevels
scmax = max(scmax,scalar_new(k,iCell))
scmin = min(scmin,scalar_new(k,iCell))
if (s_max(k,iCell) < scalar_new(k,iCell)) then
@@ -1763,8 +1781,8 @@
if (s_min(k,iCell) > scalar_new(k,iCell)) then
write(32,*) ' under - k,iCell,s_min,s_max,scalar_new ',k,iCell,s_min(k,iCell),s_max(k,iCell),scalar_new(k,iCell)
end if
- enddo
- enddo
+ end do
+ end do
write(0,*) ' scmin, scmax new out ',scmin,scmax
write(0,*) ' icell_min, k_min ',icellmax, kmax
@@ -1773,8 +1791,8 @@
! the update should be positive definite. but roundoff can sometimes leave small negative values
! hence the enforcement of PD in the copy back to the model state.
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
+ do iCell = 1, nCells
+ do k=1, nVertLevels
s_new % scalars % array(iScalar,k,iCell) = max(0.0_RKIND,scalar_new(k,iCell))
end do
end do
@@ -1814,7 +1832,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, iq
real (kind=RKIND) :: flux, workpv
- integer :: nCells, nEdges, nVertices, nVertLevels, nCellsSolve
+ integer :: nCells, nEdges, nVertices, nVertLevels, nCellsSolve, nEdgesSolve
real (kind=RKIND) :: h_mom_eddy_visc2, v_mom_eddy_visc2, h_mom_eddy_visc4
real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, h_theta_eddy_visc4
real (kind=RKIND) :: u_diffusion
@@ -1962,6 +1980,7 @@
nVertLevels = grid % nVertLevels
nVertices = grid % nVertices
nCellsSolve = grid % nCellsSolve
+ nEdgesSolve = grid % nEdgesSolve
h_mom_eddy_visc2 = config_h_mom_eddy_visc2
! h_mom_eddy_visc4 = config_h_mom_eddy_visc4
@@ -1970,6 +1989,7 @@
! h_theta_eddy_visc4 = config_h_theta_eddy_visc4
v_theta_eddy_visc2 = config_v_theta_eddy_visc2
+ nEdgesOnCell => grid % nEdgesOnCell % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
@@ -1989,7 +2009,7 @@
do iCell = 1, nCells
d_diag(:) = 0.
d_off_diag(:) = 0.
- do iEdge = 1, grid % nEdgesOnCell % array (iCell)
+ do iEdge = 1, nEdgesOnCell(iCell)
do k=1, nVertLevels
d_diag(k) = d_diag(k) + defc_a(iEdge,iCell)*u(k,EdgesOnCell(iEdge,iCell)) &
- defc_b(iEdge,iCell)*v(k,EdgesOnCell(iEdge,iCell))
@@ -2035,7 +2055,7 @@
! accumulate horizontal mass-flux
- do iEdge=1,grid % nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,nVertLevels
@@ -2067,7 +2087,7 @@
! Compute u (normal) velocity tendency for each edge (cell face)
!
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2189,7 +2209,7 @@
if (delsq_horiz_mixing) then
if ((h_mom_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1, nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -2216,7 +2236,7 @@
else if ( config_horiz_mixing == "2d_smagorinsky") then
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1, nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -2252,7 +2272,7 @@
delsq_u(:,:) = 0.0
- do iEdge=1,grid % nEdges
+ do iEdge=1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -2301,7 +2321,7 @@
end do
end do
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
@@ -2340,7 +2360,7 @@
if (config_mix_full) then
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2364,7 +2384,7 @@
else ! idealized cases where we mix on the perturbation from the initial 1-D state
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2402,7 +2422,7 @@
! add in mixing for u
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdgesSolve
do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) + tend_u_euler(k,iEdge)
end do
@@ -2422,7 +2442,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
flux = dvEdge(iEdge) * (fzm(k)*ru(k,iEdge) + fzp(k)*ru(k-1,iEdge) ) &
*(w(k,cell1) + w(k,cell2))*0.5
tend_w(k,cell1) = tend_w(k,cell1) - flux
@@ -2438,7 +2458,7 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
ru_edge_w(k) = fzm(k)*ru(k,iEdge) + fzp(k)*ru(k-1,iEdge)
end do
@@ -2448,13 +2468,13 @@
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,ru_edge_w(k))*adv_coefs_3rd(i,iEdge)
flux_arr(k) = flux_arr(k) + scalar_weight* w(k,iCell)
end do
end do
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
tend_w(k,cell1) = tend_w(k,cell1) - ru_edge_w(k)*flux_arr(k)
tend_w(k,cell2) = tend_w(k,cell2) + ru_edge_w(k)*flux_arr(k)
end do
@@ -2469,16 +2489,16 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
d2fdx2_cell1 = deriv_two(1,1,iEdge) * w(k,cell1)
d2fdx2_cell2 = deriv_two(1,2,iEdge) * w(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
- if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
+ do i=1, nEdgesOnCell(cell1)
+ if ( grid % CellsOnCell % array (i,cell1) <= nCells) &
d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * w(k,grid % CellsOnCell % array (i,cell1))
end do
- do i=1, grid % nEdgesOnCell % array (cell2)
- if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
+ do i=1, nEdgesOnCell(cell2)
+ if ( grid % CellsOnCell % array (i,cell2) <= nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * w(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -2497,7 +2517,7 @@
if (curvature) then
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
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. &
@@ -2524,14 +2544,14 @@
if ((h_mom_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed")) then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
! horizontal flux divergence of the gradient (i.e. del^2)
! note, for w, even though we use theta_* local scratch variables
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
theta_turb_flux = h_mom_eddy_visc2*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
@@ -2544,12 +2564,12 @@
else if (config_horiz_mixing == "2d_smagorinsky") then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
theta_turb_flux = 0.25*(kdiff(k,cell1)+kdiff(k,cell2)+kdiff(k-1,cell1)+kdiff(k,cell2)) &
*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
@@ -2576,10 +2596,10 @@
delsq_theta(:,:) = 0.
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=2,grid % nVertLevels
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=2,nVertLevels
delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*0.5*(rho_edge(k,iEdge)+rho_edge(k-1,iEdge))*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*0.5*(rho_edge(k,iEdge)+rho_edge(k-1,iEdge))*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
end do
@@ -2592,12 +2612,12 @@
end do
end do
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=2,grid % nVertLevels
+ do k=2,nVertLevels
theta_turb_flux = h_mom_eddy_visc4*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel4(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
@@ -2673,7 +2693,7 @@
if ( v_mom_eddy_visc2 > 0.0 ) then
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k=2,nVertLevels
tend_w_euler(k,iCell) = tend_w_euler(k,iCell) + v_mom_eddy_visc2*0.5*(rho_zz(k,iCell)+rho_zz(k-1,iCell))*( &
(w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
@@ -2688,7 +2708,7 @@
! add in mixing terms for w
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k=2,nVertLevels
tend_w(k,iCell) = tend_w(k,iCell) + tend_w_euler(k,iCell)
end do
@@ -2710,7 +2730,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
flux = dvEdge(iEdge) * ru(k,iEdge) * ( 0.5*(theta_m(k,cell1) + theta_m(k,cell2)) )
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
tend_theta(k,cell2) = tend_theta(k,cell2) + flux
@@ -2728,13 +2748,13 @@
flux_arr(:) = 0.
do i=1,nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,ru(k,iEdge))*adv_coefs_3rd(i,iEdge)
flux_arr(k) = flux_arr(k) + scalar_weight* theta_m(k,iCell)
end do
end do
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
tend_theta(k,cell1) = tend_theta(k,cell1) - ru(k,iEdge)*flux_arr(k)
tend_theta(k,cell2) = tend_theta(k,cell2) + ru(k,iEdge)*flux_arr(k)
end do
@@ -2747,17 +2767,17 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta_m(k,cell1)
d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta_m(k,cell2)
- do i=1, grid % nEdgesOnCell % array (cell1)
+ do i=1, nEdgesOnCell(cell1)
if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell1))
end do
- do i=1, grid % nEdgesOnCell % array (cell2)
+ do i=1, nEdgesOnCell(cell2)
if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta_m(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -2787,12 +2807,12 @@
if (delsq_horiz_mixing) then
if ( (h_theta_eddy_visc2 > 0.0) .and. (config_horiz_mixing == "2d_fixed") ) then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
theta_turb_flux = h_theta_eddy_visc2*prandtl_inv*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
@@ -2805,12 +2825,12 @@
else if ( ( config_horiz_mixing == "2d_smagorinsky") ) then
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl_inv &
*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel2(iEdge)
@@ -2832,10 +2852,10 @@
delsq_theta(:,:) = 0.
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevels
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*rho_edge(k,iEdge)*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*rho_edge(k,iEdge)*(theta_m(k,cell2) - theta_m(k,cell1))/dcEdge(iEdge)
end do
@@ -2848,12 +2868,12 @@
end do
end do
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,grid % nVertLevels
+ do k=1,nVertLevels
theta_turb_flux = h_theta_eddy_visc4*prandtl_inv*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
theta_turb_flux = theta_turb_flux * meshScalingDel4(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
@@ -2924,7 +2944,7 @@
if (config_mix_full) then
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k=2,nVertLevels-1
z1 = zgrid(k-1,iCell)
z2 = zgrid(k ,iCell)
@@ -2943,7 +2963,7 @@
else ! idealized cases where we mix on the perturbation from the initial 1-D state
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k=2,nVertLevels-1
z1 = zgrid(k-1,iCell)
z2 = zgrid(k ,iCell)
@@ -2966,7 +2986,7 @@
end if ! compute theta mixing on first rk_step
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, nCellsSolve
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell) + tend_theta_euler(k,iCell)
end do
@@ -2996,7 +3016,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, eoe, i
real (kind=RKIND) :: h_vertex, r
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, &
@@ -3041,10 +3061,11 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
!
! Compute height on cell edges at velocity locations
@@ -3141,7 +3162,7 @@
end do
do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
iCell = cellsOnVertex(i,iVertex)
do k = 1,nVertLevels
ke(k,iCell) = ke(k,iCell) + (1.-ke_fact)*kiteAreasOnVertex(i, iVertex) * ke_vertex(k, iVertex) / areaCell(iCell)
@@ -3173,7 +3194,7 @@
do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex = 0.0
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
@@ -3189,7 +3210,7 @@
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
iEdge = edgesOnVertex(i,iVertex)
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
@@ -3203,7 +3224,7 @@
!
pv_cell(:,:) = 0.0
do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
iCell = cellsOnVertex(i,iVertex)
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
@@ -3269,25 +3290,33 @@
type (mesh_type), intent(inout) :: grid
integer :: k,iCell,iEdge,iCell1,iCell2, cell1, cell2, coef_3rd_order
+ integer :: nCells, nEdges, nVertLevels
real (kind=RKIND) :: p0, rcv, flux
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+ cellsOnEdge => grid % cellsOnEdge % array
+
coef_3rd_order = config_coef_3rd_order
if(config_theta_adv_order /=3) coef_3rd_order = 0
rcv = rgas / (cp-rgas)
p0 = 1.e5 ! this should come from somewhere else...
- do iCell=1,grid%nCells
- do k=1,grid%nVertLevels
+ do iCell=1,nCells
+ do k=1,nVertLevels
state % theta_m % array(k,iCell) = diag % theta % array(k,iCell) * (1._RKIND + rvord * state % scalars % array(state % index_qv,k,iCell))
state % rho_zz % array(k,iCell) = diag % rho % array(k,iCell) / grid % zz % array(k,iCell)
end do
end do
- do iEdge = 1, grid % nEdges
- iCell1 = grid % cellsOnEdge % array(1,iEdge)
- iCell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevels
+ do iEdge = 1, nEdges
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
diag % ru % array(k,iEdge) = 0.5 * state % u % array(k,iEdge) * (state % rho_zz % array(k,iCell1) + state % rho_zz % array(k,iCell2))
end do
end do
@@ -3295,10 +3324,10 @@
! Compute rw (i.e. rho_zz * omega) from rho_zz, w, and ru.
! We are reversing the procedure we use in subroutine atm_recover_large_step_variables.
! first, the piece that depends on w.
- do iCell=1,grid%nCells
+ do iCell=1,nCells
diag % rw % array(1,iCell) = 0.
diag % rw % array(grid%nVertLevels+1,iCell) = 0.
- do k=2,grid%nVertLevels
+ do k=2,nVertLevels
diag % rw % array(k,iCell) = state % w % array(k,iCell) &
* (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell)) &
* (grid % fzp % array(k) * grid % zz % array(k-1,iCell) + grid % fzm % array(k) * grid % zz % array(k,iCell))
@@ -3306,10 +3335,10 @@
end do
! next, the piece that depends on ru
- do iEdge=1,grid%nEdges
- cell1 = grid % CellsOnEdge % array(1,iEdge)
- cell2 = grid % CellsOnEdge % array(2,iEdge)
- do k = 2, grid % nVertLevels
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k = 2, nVertLevels
flux = (grid % fzm % array(k) * diag % ru % array(k,iEdge)+grid % fzp % array(k) * diag % ru % array(k-1,iEdge))
diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
+ (grid % zb % array(k,2,iEdge) + coef_3rd_order * sign(1.0_RKIND,flux) * grid % zb3 % array(k,2,iEdge))*flux &
@@ -3320,33 +3349,33 @@
end do
end do
- do iCell = 1, grid % nCells
- do k=1,grid % nVertLevels
+ do iCell = 1, nCells
+ do k=1,nVertLevels
diag % rho_p % array(k,iCell) = state % rho_zz % array(k,iCell) - diag % rho_base % array(k,iCell)
end do
end do
- do iCell = 1, grid % nCells
- do k=1,grid % nVertLevels
+ do iCell = 1, nCells
+ do k=1,nVertLevels
diag % rtheta_base % array(k,iCell) = diag % theta_base % array(k,iCell) * diag % rho_base % array(k,iCell)
end do
end do
- do iCell = 1, grid % nCells
- do k=1,grid % nVertLevels
+ do iCell = 1, nCells
+ do k=1,nVertLevels
diag % rtheta_p % array(k,iCell) = state % theta_m % array(k,iCell) * diag % rho_p % array(k,iCell) &
+ diag % rho_base % array(k,iCell) * (state % theta_m % array(k,iCell) - diag % theta_base % array(k,iCell))
end do
end do
- do iCell=1,grid % nCells
- do k=1,grid % nVertLevels
+ do iCell=1,nCells
+ do k=1,nVertLevels
diag % exner % array(k,iCell) = (grid % zz % array(k,iCell) * (rgas/p0) * (diag % rtheta_p % array(k,iCell) + diag % rtheta_base % array(k,iCell)))**rcv
end do
end do
- do iCell=1,grid % nCells
- do k=1,grid % nVertLevels
+ do iCell=1,nCells
+ do k=1,nVertLevels
diag % pressure_p % array(k,iCell) = grid % zz % array(k,iCell) * rgas &
* ( diag % exner % array(k,iCell) * diag % rtheta_p % array(k,iCell) &
+ diag % rtheta_base % array(k,iCell) * (diag % exner % array(k,iCell) - diag % exner_base % array(k,iCell)) &
</font>
</pre>