<p><b>duda</b> 2010-03-03 15:07:29 -0700 (Wed, 03 Mar 2010)</p><p>Begin using the super-array functionality of the registry to<br>
manage moisture arrays. <br>
<br>
Initially add three moisture arrays, qv, qc, and qr, which <br>
will be bundled together in memory into the scalars array;<br>
these moisture arrays can be accessed individually, as <br>
<br>
scalars(index_qv,:,:)<br>
scalars(index_qc,:,:) <br>
scalars(index_qr,:,:) <br>
<br>
or they can be accessed collectively as the array <br>
<br>
scalars(moist_start:moist_end,:,:)<br>
<br>
where the constants index_qv, index_qv, index_qr, moist_start,<br>
and moist_end are automatically defined by the registry.<br>
<br>
Also, convert to the use of the term 'scalars' rather than <br>
'tracers' in the code.<br>
<br>
<br>
M src/module_test_cases.F<br>
M src/module_time_integration.F<br>
M Registry/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2010-03-03 01:16:54 UTC (rev 122)
+++ branches/hyd_model/Registry/Registry 2010-03-03 22:07:29 UTC (rev 123)
@@ -37,7 +37,7 @@
dim FIFTEEN 15
dim TWENTYONE 21
dim nVertLevels nVertLevels
-dim nTracers nTracers
+#dim nTracers nTracers
dim nVertLevelsP1 nVertLevels+1
#
@@ -102,7 +102,10 @@
var real u ( nVertLevels nEdges Time ) iro u - -
var real theta ( nVertLevels nCells Time ) iro theta - -
var real surface_pressure ( nCells Time ) iro surface_pressure - -
-var real tracers ( nTracers nVertLevels nCells Time ) iro tracers - -
+var real qv ( nVertLevels nCells Time ) iro qv scalars moist
+var real qc ( nVertLevels nCells Time ) iro qc scalars moist
+var real qr ( nVertLevels nCells Time ) iro qr scalars moist
+#var real tracers ( nTracers nVertLevels nCells Time ) iro tracers - -
# state variables diagnosed from prognostic state
var real h ( nVertLevels nCells Time ) ro h - -
@@ -140,7 +143,10 @@
var real h_edge_old ( nVertLevels nEdges ) - h_edge_old - -
var real h_old ( nVertLevels nCells ) - h_old - -
var real pressure_old ( nVertLevelsP1 nCells ) - pressure_old - -
-var real tracers_old ( nTracers nVertLevels nCells ) - tracers_old - -
+var real qv_old ( nVertLevels nCells ) - qv_old scalars_old moist_old
+var real qc_old ( nVertLevels nCells ) - qc_old scalars_old moist_old
+var real qr_old ( nVertLevels nCells ) - qr_old scalars_old moist_old
+#var real tracers_old ( nTracers nVertLevels nCells ) - tracers_old - -
# Space needed for advection
var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
Modified: branches/hyd_model/src/module_test_cases.F
===================================================================
--- branches/hyd_model/src/module_test_cases.F 2010-03-03 01:16:54 UTC (rev 122)
+++ branches/hyd_model/src/module_test_cases.F 2010-03-03 22:07:29 UTC (rev 123)
@@ -139,7 +139,7 @@
alpha => state % alpha % array
geopotential => state % geopotential % array
h => state % h % array
- scalars => state % tracers % array
+ scalars => state % scalars % array
scalars = 0.
@@ -341,22 +341,24 @@
theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
end do
+! When initializing a scalar, be sure not to put unreasonably large values
+! into indices in the moist class
! scalars(2,:,:) = 1. ! transport test
! scalars(2,:,:) = theta ! transport test
- if (grid % nTracers >= 2) then
- scalars(2,:,:) = 0.0
- do iCell=1,grid%nCells
- r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
- if (r < a/3.0) then
- do k=1,grid%nVertLevels
- scalars(2,k,iCell) = (1.0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
- end do
- end if
- end do
- end if
- if (grid % nTracers >= 3) scalars(3,:,:) = theta + 100. ! transport test
- if (grid % nTracers >= 4) scalars(4,:,:) = theta + 200. ! transport test
- if (grid % nTracers >= 5) scalars(5,:,:) = theta + 300. ! transport test
+! if (num_scalars >= 2) then
+! scalars(2,:,:) = 0.0
+! do iCell=1,grid%nCells
+! r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
+! if (r < a/3.0) then
+! do k=1,grid%nVertLevels
+! scalars(2,k,iCell) = (1.0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
+! end do
+! end if
+! end do
+! end if
+! if (num_scalars >= 3) scalars(3,:,:) = theta + 100. ! transport test
+! if (num_scalars >= 4) scalars(4,:,:) = theta + 200. ! transport test
+! if (num_scalars >= 5) scalars(5,:,:) = theta + 300. ! transport test
end subroutine hyd_test_case_1
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2010-03-03 01:16:54 UTC (rev 122)
+++ branches/hyd_model/src/module_time_integration.F 2010-03-03 22:07:29 UTC (rev 123)
@@ -73,8 +73,8 @@
logical, parameter :: debug = .false.
logical, parameter :: debug_mass_conservation = .true.
- real (kind=RKIND) :: domain_mass, tracer_mass, tracer_min, tracer_max
- real (kind=RKIND) :: global_domain_mass, global_tracer_mass, global_tracer_min, global_tracer_max
+ real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
+ real (kind=RKIND) :: global_domain_mass, global_scalar_mass, global_scalar_min, global_scalar_max
!
! Initialize time_levs(2) with state at current time
@@ -176,13 +176,13 @@
block => domain % blocklist
do while (associated(block))
!
- ! Note: Tracers in new time level shouldn't be overwritten, since their provisional values
+ ! Note: Scalars in new time level shouldn't be overwritten, since their provisional values
! from the previous RK step are needed to compute new scalar tendencies in advance_scalars.
- ! A cleaner way of preserving tracers should be added in future.
+ ! A cleaner way of preserving scalars should be added in future.
!
- block % mesh % tracers_old % array(:,:,:) = block % time_levs(2) % state % tracers % array(:,:,:)
+ block % mesh % scalars_old % array(:,:,:) = block % time_levs(2) % state % scalars % array(:,:,:)
call copy_state( block % time_levs(1) % state, block % time_levs(2) % state )
- block % time_levs(2) % state % tracers % array(:,:,:) = block % mesh % tracers_old % array(:,:,:)
+ block % time_levs(2) % state % scalars % array(:,:,:) = block % mesh % scalars_old % array(:,:,:)
block => block % next
end do
@@ -282,11 +282,11 @@
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -318,10 +318,10 @@
if(debug .or. debug_mass_conservation) then
domain_mass = 0.
- tracer_mass = 0.
+ scalar_mass = 0.
block => domain % blocklist
- tracer_min = block % time_levs(2) % state % tracers % array (2,1,1)
- tracer_max = block % time_levs(2) % state % tracers % array (2,1,1)
+ scalar_min = block % time_levs(2) % state % scalars % array (2,1,1)
+ scalar_max = block % time_levs(2) % state % scalars % array (2,1,1)
do while(associated(block))
do iCell = 1, block % mesh % nCellsSolve
domain_mass = domain_mass + block % time_levs(2) % state % surface_pressure % array (iCell) * &
@@ -329,23 +329,23 @@
- block % time_levs(2) % state % pressure % array (block % mesh % nVertLevels + 1, 1) * &
block % mesh % areaCell % array (iCell)
do k=1, block % mesh % nVertLevelsSolve
- tracer_mass = tracer_mass - block % time_levs(2) % state % tracers % array (2,k,iCell) * &
+ scalar_mass = scalar_mass - block % time_levs(2) % state % scalars % array (2,k,iCell) * &
block % time_levs(2) % state % h % array (k,iCell) * &
block % mesh % dnw % array (k) * &
block % mesh % areaCell % array (iCell)
- tracer_min = min(tracer_min,block % time_levs(2) % state % tracers % array (2,k,iCell))
- tracer_max = max(tracer_max,block % time_levs(2) % state % tracers % array (2,k,iCell))
+ scalar_min = min(scalar_min,block % time_levs(2) % state % scalars % array (2,k,iCell))
+ scalar_max = max(scalar_max,block % time_levs(2) % state % scalars % array (2,k,iCell))
end do
end do
block => block % next
end do
call dmpar_sum_real(domain % dminfo, domain_mass, global_domain_mass)
- call dmpar_sum_real(domain % dminfo, tracer_mass, global_tracer_mass)
- call dmpar_min_real(domain % dminfo, tracer_min, global_tracer_min)
- call dmpar_max_real(domain % dminfo, tracer_max, global_tracer_max)
+ call dmpar_sum_real(domain % dminfo, scalar_mass, global_scalar_mass)
+ call dmpar_min_real(domain % dminfo, scalar_min, global_scalar_min)
+ call dmpar_max_real(domain % dminfo, scalar_max, global_scalar_max)
write(0,*) ' mass in the domain = ',global_domain_mass
- write(0,*) ' tracer mass in the domain = ',global_tracer_mass
- write(0,*) ' tracer_min, tracer_max ',global_tracer_min, global_tracer_max
+ write(0,*) ' scalar mass in the domain = ',global_scalar_mass
+ write(0,*) ' scalar_min, scalar_max ',global_scalar_min, global_scalar_max
end if
@@ -372,23 +372,21 @@
integer :: iEdge, iCell, k, cell1, cell2, iq
- integer :: nCells, nEdges, nVertLevels, nTracers
+ integer :: nCells, nEdges, nVertLevels
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
- nTracers = grid % nTracers
grid % qtot % array = 0.
grid % cqu % array = 1.
- if(nTracers > 0) then ! we assume these are all moisture variables at present
+ if (num_scalars > 0) then
do iCell = 1, nCells
do k = 1, nVertLevels
-! do iq = 1, nTracers
- do iq = 1, 1
- grid % qtot % array(k,iCell) = grid % qtot % array(k,iCell) + s % tracers % array (iq, k, iCell)
+ do iq = moist_start, moist_end
+ grid % qtot % array(k,iCell) = grid % qtot % array(k,iCell) + s % scalars % array (iq, k, iCell)
end do
end do
end do
@@ -432,7 +430,7 @@
real (kind=RKIND), dimension(:), pointer :: surface_pressure, dbn, dnu, dnw
integer :: iEdge, iCell, k, cell1, cell2, iq
- integer :: nCells, nEdges, nVertLevels, nTracers
+ integer :: nCells, nEdges, nVertLevels
real (kind=RKIND) :: p0,tm,ptop,ptmp
@@ -443,7 +441,7 @@
surface_pressure => s % surface_pressure % array
alpha => s % alpha % array
geopotential => s % geopotential % array
- scalar => s % tracers % array
+ scalar => s % scalars % array
theta_old => grid % theta_old % array
u_old => grid % u_old % array
ww_old => grid % ww_old % array
@@ -460,7 +458,6 @@
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
- nTracers = grid % nTracers
@@ -1115,7 +1112,7 @@
qtot => grid % qtot % array
cqu => grid % cqu % array
ww => s % ww % array
- scalar => s % tracers % array
+ scalar => s % scalars % array
dnw => grid % dnw % array
dbn => grid % dbn % array
@@ -1312,8 +1309,8 @@
type (grid_meta), intent(in) :: grid
real (kind=RKIND) :: dt
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge, d2fdx2_cell1, d2fdx2_cell2
+ integer :: i, iCell, iEdge, k, iScalar, cell1, cell2
+ real (kind=RKIND) :: flux, scalar_edge, d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
@@ -1321,7 +1318,7 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( grid % nTracers, grid % nVertLevels + 1 ) :: wdtn
+ real (kind=RKIND), dimension( num_scalars, grid % nVertLevels + 1 ) :: wdtn
integer :: nVertLevels
real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw
@@ -1331,14 +1328,14 @@
if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
if (config_scalar_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
- scalar_old => s_old % tracers % array
- scalar_new => s_new % tracers % array
+ scalar_old => s_old % scalars % array
+ scalar_new => s_new % scalars % array
deriv_two => grid % deriv_two % array
uhAvg => grid % uhAvg % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
cellsOnEdge => grid % cellsOnEdge % array
- scalar_tend => tend % tracers % array
+ scalar_tend => tend % scalars % array
h_old => s_old % h % array
h_new => s_new % h % array
wwAvg => grid % wwAvg % array
@@ -1365,11 +1362,11 @@
cell2 = cellsOnEdge(2,iEdge)
if (cell1 > 0 .and. cell2 > 0) then
do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2))
- flux = uhAvg(k,iEdge) * dvEdge(iEdge) * tracer_edge
- scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ do iScalar=1,num_scalars
+ scalar_edge = 0.5 * (scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2))
+ flux = uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_edge
+ scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) - flux/areaCell(cell1)
+ scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) + flux/areaCell(cell2)
end do
end do
end if
@@ -1384,28 +1381,28 @@
do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ do iScalar=1,num_scalars
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ deriv_two(i+1,1,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
if (uhAvg(k,iEdge) > 0) then
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
else
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
@@ -1413,16 +1410,16 @@
! old version of the above code, with coef_3rd_order assumed to be 1.0
! if (uhAvg(k,iEdge) > 0) then
! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
! else
! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
-! 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+! 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
! end if
- scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) - flux/areaCell(cell1)
+ scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) + flux/areaCell(cell2)
end do
end do
@@ -1438,26 +1435,26 @@
do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ do iScalar=1,num_scalars
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ deriv_two(i+1,1,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
- scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) - flux/areaCell(cell1)
+ scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) + flux/areaCell(cell2)
end do
end do
end if
@@ -1474,16 +1471,16 @@
wdtn(:,1) = 0.
do k = 2, nVertLevels
- do iTracer=1,grid % nTracers
- wdtn(iTracer,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iTracer,k,iCell)+fnp(k)*scalar_new(iTracer,k-1,iCell))
+ do iScalar=1,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
wdtn(:,nVertLevels+1) = 0.
do k=1,grid % nVertLevelsSolve
- do iTracer=1,grid % nTracers
- scalar_new(iTracer,k,iCell) = ( scalar_old(iTracer,k,iCell)*h_old(k,iCell) &
- + dt*( scalar_tend(iTracer,k,iCell) -rdnw(k)*(wdtn(iTracer,k+1)-wdtn(iTracer,k)) ) )/h_new(k,iCell)
+ do iScalar=1,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)
end do
end do
@@ -1512,8 +1509,8 @@
type (dm_info), intent(in) :: dminfo
type (exchange_list), pointer :: cellsToSend, cellsToRecv
- integer :: i, iCell, iEdge, k, iTracer, cell_upwind, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge, d2fdx2_cell1, d2fdx2_cell2
+ integer :: i, iCell, iEdge, k, iScalar, cell_upwind, cell1, cell2
+ real (kind=RKIND) :: flux, scalar_edge, d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND) :: fdir, flux_upwind, h_flux_upwind, s_upwind
real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
@@ -1522,10 +1519,10 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( grid % nTracers, grid % nEdges) :: h_flux
- real (kind=RKIND), dimension( grid % nTracers, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
- real (kind=RKIND), dimension( grid % nTracers, grid % nCells, 2 ) :: scale_out, scale_in
- real (kind=RKIND), dimension( grid % nTracers ) :: s_max, s_min, s_max_update, s_min_update
+ real (kind=RKIND), dimension( num_scalars, grid % nEdges) :: h_flux
+ real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( num_scalars ) :: s_max, s_min, s_max_update, s_min_update
integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
@@ -1533,14 +1530,14 @@
real (kind=RKIND), parameter :: eps=1.e-20
real (kind=RKIND) :: coef_3rd_order
- scalar_old => s_old % tracers % array
- scalar_new => s_new % tracers % array
+ scalar_old => s_old % scalars % array
+ scalar_new => s_new % scalars % array
deriv_two => grid % deriv_two % array
uhAvg => grid % uhAvg % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
cellsOnEdge => grid % cellsOnEdge % array
- scalar_tend => tend % tracers % array
+ scalar_tend => tend % scalars % array
h_old => s_old % h % array
h_new => s_new % h % array
wwAvg => grid % wwAvg % array
@@ -1580,21 +1577,21 @@
if (k < grid % nVertLevels) then
cell_upwind = k
if (wwAvg(k+1,iCell) >= 0) cell_upwind = k+1
- do iTracer=1,grid % nTracers
- v_flux(iTracer,iCell,km0) = dt * wwAvg(k+1,iCell) * &
- (fnm(k+1) * scalar_new(iTracer,k+1,iCell) + fnp(k+1) * scalar_new(iTracer,k,iCell))
- v_flux_upwind(iTracer,iCell,km0) = dt * wwAvg(k+1,iCell) * scalar_old(iTracer,cell_upwind,iCell)
- v_flux(iTracer,iCell,km0) = v_flux(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km0)
-! v_flux(iTracer,iCell,km0) = 0. ! use only upwind - for testing
- s_update(iTracer,iCell,km0) = scalar_old(iTracer,k,iCell) * h_old(k,iCell) &
- - rdnw(k) * (v_flux_upwind(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km1))
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * &
+ (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
+ v_flux_upwind(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) * scalar_old(iScalar,cell_upwind,iCell)
+ v_flux(iScalar,iCell,km0) = v_flux(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km0)
+! v_flux(iScalar,iCell,km0) = 0. ! use only upwind - for testing
+ s_update(iScalar,iCell,km0) = scalar_old(iScalar,k,iCell) * h_old(k,iCell) &
+ - rdnw(k) * (v_flux_upwind(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km1))
end do
else
- do iTracer=1,grid % nTracers
- v_flux(iTracer,iCell,km0) = 0.
- v_flux_upwind(iTracer,iCell,km0) = 0.
- s_update(iTracer,iCell,km0) = scalar_old(iTracer,k,iCell) * h_old(k,iCell) &
- - rdnw(k) * (v_flux_upwind(iTracer,iCell,km0) - v_flux_upwind(iTracer,iCell,km1))
+ do iScalar=1,num_scalars
+ v_flux(iScalar,iCell,km0) = 0.
+ v_flux_upwind(iScalar,iCell,km0) = 0.
+ s_update(iScalar,iCell,km0) = scalar_old(iScalar,k,iCell) * h_old(k,iCell) &
+ - rdnw(k) * (v_flux_upwind(iScalar,iCell,km0) - v_flux_upwind(iScalar,iCell,km1))
end do
end if
@@ -1610,14 +1607,14 @@
if (cell1 > 0 .and. cell2 > 0) then
cell_upwind = cell2
if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2))
- h_flux(iTracer,iEdge) = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * tracer_edge
- h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iTracer,k,cell_upwind)
- h_flux(iTracer,iEdge) = h_flux(iTracer,iEdge) - h_flux_upwind
-! h_flux(iTracer,iEdge) = 0. ! use only upwind - for testing
- s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
- s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
+ do iScalar=1,num_scalars
+ scalar_edge = 0.5 * (scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2))
+ h_flux(iScalar,iEdge) = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_edge
+ h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iScalar,k,cell_upwind)
+ h_flux(iScalar,iEdge) = h_flux(iScalar,iEdge) - h_flux_upwind
+! h_flux(iScalar,iEdge) = 0. ! use only upwind - for testing
+ s_update(iScalar,cell1,km0) = s_update(iScalar,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
+ s_update(iScalar,cell2,km0) = s_update(iScalar,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
end do
end if
end do
@@ -1630,39 +1627,39 @@
if (cell1 > 0 .and. cell2 > 0) then
cell_upwind = cell2
if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
- do iTracer=1,grid % nTracers
+ do iScalar=1,num_scalars
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ deriv_two(i+1,1,iEdge) * scalar_new(iScalar,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) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
if (uhAvg(k,iEdge) > 0) then
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
-(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
else
flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
- h_flux(iTracer,iEdge) = dt * flux
- h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iTracer,k,cell_upwind)
- h_flux(iTracer,iEdge) = h_flux(iTracer,iEdge) - h_flux_upwind
-! h_flux(iTracer,iEdge) = 0. ! use only upwind - for testing
- s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
- s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
+ h_flux(iScalar,iEdge) = dt * flux
+ h_flux_upwind = dt * uhAvg(k,iEdge) * dvEdge(iEdge) * scalar_old(iScalar,k,cell_upwind)
+ h_flux(iScalar,iEdge) = h_flux(iScalar,iEdge) - h_flux_upwind
+! h_flux(iScalar,iEdge) = 0. ! use only upwind - for testing
+ s_update(iScalar,cell1,km0) = s_update(iScalar,cell1,km0) - h_flux_upwind / grid % areaCell % array(cell1)
+ s_update(iScalar,cell2,km0) = s_update(iScalar,cell2,km0) + h_flux_upwind / grid % areaCell % array(cell2)
end do
end if
end do
@@ -1674,31 +1671,31 @@
!*************************************************************************************************************
!--- limiter - we limit horizontal and vertical fluxes on level k
-!--- (these are h fluxes contributing to level k scalars, and v flux contributing to level k, k-1 tracers)
+!--- (these are h fluxes contributing to level k scalars, and v flux contributing to level k, k-1 scalars)
do iCell=1,grid % nCells
- do iTracer=1,grid % nTracers
+ do iScalar=1,num_scalars
- s_max(iTracer) = max(scalar_old(iTracer,k,iCell), scalar_old(iTracer,kcp1,iCell), scalar_old(iTracer,kcm1,iCell))
- s_min(iTracer) = min(scalar_old(iTracer,k,iCell), scalar_old(iTracer,kcp1,iCell), scalar_old(iTracer,kcm1,iCell))
- s_max_update(iTracer) = s_update(iTracer,iCell,km0)
- s_min_update(iTracer) = s_update(iTracer,iCell,km0)
+ s_max(iScalar) = max(scalar_old(iScalar,k,iCell), scalar_old(iScalar,kcp1,iCell), scalar_old(iScalar,kcm1,iCell))
+ s_min(iScalar) = min(scalar_old(iScalar,k,iCell), scalar_old(iScalar,kcp1,iCell), scalar_old(iScalar,kcm1,iCell))
+ s_max_update(iScalar) = s_update(iScalar,iCell,km0)
+ s_min_update(iScalar) = s_update(iScalar,iCell,km0)
! add in vertical flux to get max and min estimate
- s_max_update(iTracer) = s_max_update(iTracer) &
- - rdnw(k) * (max(0.,v_flux(iTracer,iCell,km0)) - min(0.,v_flux(iTracer,iCell,km1)))
- s_min_update(iTracer) = s_min_update(iTracer) &
- - rdnw(k) * (min(0.,v_flux(iTracer,iCell,km0)) - max(0.,v_flux(iTracer,iCell,km1)))
+ s_max_update(iScalar) = s_max_update(iScalar) &
+ - rdnw(k) * (max(0.,v_flux(iScalar,iCell,km0)) - min(0.,v_flux(iScalar,iCell,km1)))
+ s_min_update(iScalar) = s_min_update(iScalar) &
+ - rdnw(k) * (min(0.,v_flux(iScalar,iCell,km0)) - max(0.,v_flux(iScalar,iCell,km1)))
end do
do i = 1, grid % nEdgesOnCell % array(iCell) ! go around the edges of each cell
if (grid % cellsOnCell % array(i,iCell) > 0) then
- do iTracer=1,grid % nTracers
+ do iScalar=1,num_scalars
- s_max(iTracer) = max(scalar_old(iTracer,k,grid % cellsOnCell % array(i,iCell)), s_max(iTracer))
- s_min(iTracer) = min(scalar_old(iTracer,k,grid % cellsOnCell % array(i,iCell)), s_min(iTracer))
+ s_max(iScalar) = max(scalar_old(iScalar,k,grid % cellsOnCell % array(i,iCell)), s_max(iScalar))
+ s_min(iScalar) = min(scalar_old(iScalar,k,grid % cellsOnCell % array(i,iCell)), s_min(iScalar))
iEdge = grid % EdgesOnCell % array (i,iCell)
if (iCell == cellsOnEdge(1,iEdge)) then
@@ -1706,9 +1703,9 @@
else
fdir = -1.0
end if
- flux = -fdir * h_flux(iTracer,iEdge)/grid % areaCell % array(iCell)
- s_max_update(iTracer) = s_max_update(iTracer) + max(0.,flux)
- s_min_update(iTracer) = s_min_update(iTracer) + min(0.,flux)
+ flux = -fdir * h_flux(iScalar,iEdge)/grid % areaCell % array(iCell)
+ s_max_update(iScalar) = s_max_update(iScalar) + max(0.,flux)
+ s_min_update(iScalar) = s_min_update(iScalar) + min(0.,flux)
end do
end if
@@ -1717,32 +1714,32 @@
if( config_positive_definite ) s_min(:) = 0.
- do iTracer=1,grid % nTracers
- scale_out (iTracer,iCell,km0) = 1.
- scale_in (iTracer,iCell,km0) = 1.
- s_max_update (iTracer) = s_max_update (iTracer) / h_new (k,iCell)
- s_min_update (iTracer) = s_min_update (iTracer) / h_new (k,iCell)
- s_upwind = s_update(iTracer,iCell,km0) / h_new(k,iCell)
- if ( s_max_update(iTracer) > s_max(iTracer) .and. config_monotonic) &
- scale_in (iTracer,iCell,km0) = max(0.,(s_max(iTracer)-s_upwind)/(s_max_update(iTracer)-s_upwind+eps))
- if ( s_min_update(iTracer) < s_min(iTracer) ) &
- scale_out (iTracer,iCell,km0) = max(0.,(s_upwind-s_min(iTracer))/(s_upwind-s_min_update(iTracer)+eps))
+ do iScalar=1,num_scalars
+ scale_out (iScalar,iCell,km0) = 1.
+ scale_in (iScalar,iCell,km0) = 1.
+ s_max_update (iScalar) = s_max_update (iScalar) / h_new (k,iCell)
+ s_min_update (iScalar) = s_min_update (iScalar) / h_new (k,iCell)
+ s_upwind = s_update(iScalar,iCell,km0) / h_new(k,iCell)
+ if ( s_max_update(iScalar) > s_max(iScalar) .and. config_monotonic) &
+ scale_in (iScalar,iCell,km0) = max(0.,(s_max(iScalar)-s_upwind)/(s_max_update(iScalar)-s_upwind+eps))
+ if ( s_min_update(iScalar) < s_min(iScalar) ) &
+ scale_out (iScalar,iCell,km0) = max(0.,(s_upwind-s_min(iScalar))/(s_upwind-s_min_update(iScalar)+eps))
end do
end do ! end loop over cells to compute scale factor
call dmpar_exch_halo_field2dReal(dminfo, scale_out(:,:,1), &
- grid % nTracers, grid % nCells, &
+ num_scalars, grid % nCells, &
cellsToSend, cellsToRecv)
call dmpar_exch_halo_field2dReal(dminfo, scale_out(:,:,2), &
- grid % nTracers, grid % nCells, &
+ num_scalars, grid % nCells, &
cellsToSend, cellsToRecv)
call dmpar_exch_halo_field2dReal(dminfo, scale_in(:,:,1), &
- grid % nTracers, grid % nCells, &
+ num_scalars, grid % nCells, &
cellsToSend, cellsToRecv)
call dmpar_exch_halo_field2dReal(dminfo, scale_in(:,:,2), &
- grid % nTracers, grid % nCells, &
+ num_scalars, grid % nCells, &
cellsToSend, cellsToRecv)
! rescale the horizontal fluxes
@@ -1751,14 +1748,14 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
if (cell1 > 0 .and. cell2 > 0) then
- do iTracer=1,grid % nTracers
- flux = h_flux(iTracer,iEdge)
+ do iScalar=1,num_scalars
+ flux = h_flux(iScalar,iEdge)
if (flux > 0) then
- flux = flux * min(scale_out(iTracer,cell1,km0), scale_in(iTracer,cell2,km0))
+ flux = flux * min(scale_out(iScalar,cell1,km0), scale_in(iScalar,cell2,km0))
else
- flux = flux * min(scale_in(iTracer,cell1,km0), scale_out(iTracer,cell2,km0))
+ flux = flux * min(scale_in(iScalar,cell1,km0), scale_out(iScalar,cell2,km0))
end if
- h_flux(iTracer,iEdge) = flux
+ h_flux(iScalar,iEdge) = flux
end do
end if
end do
@@ -1766,14 +1763,14 @@
! rescale the vertical flux
do iCell=1,grid % nCells
- do iTracer=1,grid % nTracers
- flux = v_flux(iTracer,iCell,km1)
+ do iScalar=1,num_scalars
+ flux = v_flux(iScalar,iCell,km1)
if (flux > 0) then
- flux = flux * min(scale_out(iTracer,iCell,km0), scale_in(iTracer,iCell,km1))
+ flux = flux * min(scale_out(iScalar,iCell,km0), scale_in(iScalar,iCell,km1))
else
- flux = flux * min(scale_in(iTracer,iCell,km0), scale_out(iTracer,iCell,km1))
+ flux = flux * min(scale_in(iScalar,iCell,km0), scale_out(iScalar,iCell,km1))
end if
- v_flux(iTracer,iCell,km1) = flux
+ v_flux(iScalar,iCell,km1) = flux
end do
end do
@@ -1786,9 +1783,9 @@
do iCell=1,grid % nCells
! add in upper vertical flux that was just renormalized
- do iTracer=1,grid % nTracers
- s_update(iTracer,iCell,km0) = s_update(iTracer,iCell,km0) + rdnw(k) * v_flux(iTracer,iCell,km1)
- if (k > 1) s_update(iTracer,iCell,km1) = s_update(iTracer,iCell,km1) - rdnw(k-1)*v_flux(iTracer,iCell,km1)
+ do iScalar=1,num_scalars
+ s_update(iScalar,iCell,km0) = s_update(iScalar,iCell,km0) + rdnw(k) * v_flux(iScalar,iCell,km1)
+ if (k > 1) s_update(iScalar,iCell,km1) = s_update(iScalar,iCell,km1) - rdnw(k-1)*v_flux(iScalar,iCell,km1)
end do
end do
@@ -1796,11 +1793,11 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 > 0 .and. cell2 > 0) then
- do iTracer=1,grid % nTracers
- s_update(iTracer,cell1,km0) = s_update(iTracer,cell1,km0) - &
- h_flux(iTracer,iEdge) / grid % areaCell % array(cell1)
- s_update(iTracer,cell2,km0) = s_update(iTracer,cell2,km0) + &
- h_flux(iTracer,iEdge) / grid % areaCell % array(cell2)
+ do iScalar=1,num_scalars
+ s_update(iScalar,cell1,km0) = s_update(iScalar,cell1,km0) - &
+ h_flux(iScalar,iEdge) / grid % areaCell % array(cell1)
+ s_update(iScalar,cell2,km0) = s_update(iScalar,cell2,km0) + &
+ h_flux(iScalar,iEdge) / grid % areaCell % array(cell2)
end do
end if
end do
@@ -1808,14 +1805,14 @@
! decouple from mass
if (k > 1) then
do iCell=1,grid % nCells
- do iTracer=1,grid % nTracers
- s_update(iTracer,iCell,km1) = s_update(iTracer,iCell,km1) / h_new(k-1,iCell)
+ do iScalar=1,num_scalars
+ s_update(iScalar,iCell,km1) = s_update(iScalar,iCell,km1) / h_new(k-1,iCell)
end do
end do
do iCell=1,grid % nCells
- do iTracer=1,grid % nTracers
- scalar_new(iTracer,k-1,iCell) = s_update(iTracer,iCell,km1)
+ do iScalar=1,num_scalars
+ scalar_new(iScalar,k-1,iCell) = s_update(iScalar,iCell,km1)
end do
end do
end if
@@ -1827,8 +1824,8 @@
end do
do iCell=1,grid % nCells
- do iTracer=1,grid % nTracers
- scalar_new(iTracer,grid % nVertLevels,iCell) = s_update(iTracer,iCell,km1) / h_new(grid%nVertLevels,iCell)
+ do iScalar=1,num_scalars
+ scalar_new(iScalar,grid % nVertLevels,iCell) = s_update(iScalar,iCell,km1) / h_new(grid%nVertLevels,iCell)
end do
end do
</font>
</pre>