<p><b>mpetersen@lanl.gov</b> 2010-06-01 14:30:27 -0600 (Tue, 01 Jun 2010)</p><p>Merged trunk to z-level branch.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/Makefile
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/Makefile 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/Makefile 2010-06-01 20:30:27 UTC (rev 320)
@@ -1,7 +1,7 @@
#MODEL_FORMULATION = -DNCAR_FORMULATION
MODEL_FORMULATION = -DLANL_FORMULATION
-#EXPAND_LEVELS = -DEXPAND_LEVELS=26
+# EXPAND_LEVELS = -DEXPAND_LEVELS=26
FILE_OFFSET = -DOFFSET64BIT
#########################
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/Registry 2010-06-01 20:30:27 UTC (rev 320)
@@ -112,6 +112,7 @@
# state variables diagnosed from prognostic state
var real h ( nVertLevels nCells Time ) ro h - -
var real ww ( nVertLevelsP1 nCells Time ) ro ww - -
+var real w ( nVertLevelsP1 nCells Time ) ro w - -
var real pressure ( nVertLevelsP1 nCells Time ) ro pressure - -
var real geopotential ( nVertLevelsP1 nCells Time ) ro geopotential - -
var real alpha ( nVertLevels nCells Time ) iro alpha - -
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_advection.F 2010-06-01 20:30:27 UTC (rev 320)
@@ -25,6 +25,7 @@
! local variables
real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ real (kind=RKIND), dimension(grid % nEdges) :: xe, ye
real (kind=RKIND), dimension(grid % nCells) :: theta_abs
real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
@@ -102,50 +103,62 @@
if ( .not. do_the_cell ) cycle
+
! compute poynomial fit for this cell if all needed neighbors exist
+ if ( grid % on_a_sphere ) then
- 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
- end do
+ 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
+ end do
- theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
- xc(2), yc(2), zc(2), &
- 0., 0., 1. )
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
! angles from cell center to neighbor centers (thetav)
- do i=1,n-1
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) ip2 = 2
+
+ thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1), &
+ xc(ip2), yc(ip2), zc(ip2) )
- ip2 = i+2
- if (ip2 > n) ip2 = 2
-
- thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
- 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) )
+ end do
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
- end do
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ end do
- length_scale = 1.
- do i=1,n-1
- dl_sphere(i) = dl_sphere(i)/length_scale
- end do
+! thetat(1) = 0. ! this defines the x direction, cell center 1 ->
+ thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
+ do i=2,n-1
+ thetat(i) = thetat(i-1) + thetav(i-1)
+ end do
+
+ do i=1,n-1
+ xp(i) = cos(thetat(i)) * dl_sphere(i)
+ yp(i) = sin(thetat(i)) * dl_sphere(i)
+ end do
-! thetat(1) = 0. ! this defines the x direction, cell center 1 ->
- thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
- do i=2,n-1
- thetat(i) = thetat(i-1) + thetav(i-1)
- end do
+ else ! On an x-y plane
- do i=1,n-1
- xp(i) = cos(thetat(i)) * dl_sphere(i)
- yp(i) = sin(thetat(i)) * dl_sphere(i)
- end do
+ do i=1,n-1
+ xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+ yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+ end do
+ end if
+
+
ma = n-1
mw = grid % nEdgesOnCell % array (iCell)
@@ -244,20 +257,25 @@
yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- call arc_bisect( xv1, yv1, zv1, &
- xv2, yv2, zv2, &
- xec, yec, zec )
+ if ( grid % on_a_sphere ) then
+ call arc_bisect( xv1, yv1, zv1, &
+ xv2, yv2, zv2, &
+ xec, yec, zec )
- thetae_tmp = sphere_angle( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1), &
- xec, yec, zec )
- thetae_tmp = thetae_tmp + thetat(i)
-
- if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ thetae_tmp = sphere_angle( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1), &
+ xec, yec, zec )
+ thetae_tmp = thetae_tmp + thetat(i)
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ else
+ thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ end if
else
- thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
+ ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
end if
+
end do
! fill second derivative stencil for rk advection
@@ -266,31 +284,40 @@
iEdge = grid % EdgesOnCell % array (i,iCell)
- if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ if ( grid % on_a_sphere ) then
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- cos2t = cos(thetae(1,grid % EdgesOnCell % array (i,iCell)))
- sin2t = sin(thetae(1,grid % EdgesOnCell % array (i,iCell)))
- costsint = cos2t*sin2t
- cos2t = cos2t**2
- sin2t = sin2t**2
+ cos2t = cos(thetae(1,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(thetae(1,grid % EdgesOnCell % array (i,iCell)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
- do j=1,n
- deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
- + 2.*costsint*bmatrix(5,j) &
- + 2.*sin2t*bmatrix(6,j)
- end do
+ do j=1,n
+ deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ else
+
+ cos2t = cos(thetae(2,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(thetae(2,grid % EdgesOnCell % array (i,iCell)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+ do j=1,n
+ deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ end if
else
-
- cos2t = cos(thetae(2,grid % EdgesOnCell % array (i,iCell)))
- sin2t = sin(thetae(2,grid % EdgesOnCell % array (i,iCell)))
- costsint = cos2t*sin2t
- cos2t = cos2t**2
- sin2t = sin2t**2
-
do j=1,n
- deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
- + 2.*costsint*bmatrix(5,j) &
- + 2.*sin2t*bmatrix(6,j)
+ deriv_two(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) &
+ + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) &
+ + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
+ deriv_two(j,2,iEdge) = deriv_two(j,1,iEdge)
end do
end if
end do
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_hyd_atmos/module_time_integration.F 2010-06-01 20:30:27 UTC (rev 320)
@@ -314,6 +314,14 @@
! END RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! compute vertical velocity diagnostic
+
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_w( block % time_levs(2) % state, block % time_levs(1) % state, block % mesh, dt )
+ block => block % next
+ end do
+
if(debug) write(0,*) ' rk step complete - mass diagnostics '
if(debug .or. debug_mass_conservation) then
@@ -395,9 +403,7 @@
do k = 1, nVertLevels
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- grid % cqu % array(k,iEdge) = 1./( 1. + 0.5*(grid % qtot % array(k,cell1)+grid % qtot % array(k,cell2)) )
- end if
+ grid % cqu % array(k,iEdge) = 1./( 1. + 0.5*(grid % qtot % array(k,cell1)+grid % qtot % array(k,cell2)) )
end do
end do
@@ -684,33 +690,25 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
+ do k=1,nVertLevels
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for h_mom_eddy_visc4 == constant
- !
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
-
- delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
- end do
- end if
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
+ ! only valid for h_mom_eddy_visc4 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
+ end do
end do
delsq_circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
- do k=1,nVertLevels
- delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
- end do
- end if
- if (verticesOnEdge(2,iEdge) > 0) then
- do k=1,nVertLevels
- delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
- end do
- end if
+ do k=1,nVertLevels
+ delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
end do
do iVertex=1,nVertices
r = 1.0 / areaTriangle(iVertex)
@@ -723,16 +721,10 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
- do k=1,nVertLevels
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) + delsq_u(k,iEdge)*dvEdge(iEdge)
- end do
- end if
- if(cell2 > 0) then
- do k=1,nVertLevels
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) - delsq_u(k,iEdge)*dvEdge(iEdge)
- end do
- end if
+ do k=1,nVertLevels
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) + delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
@@ -833,16 +825,14 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) + flux
- tend_theta(k,cell2) = tend_theta(k,cell2) - flux
- end do
+ do k=1,grid % nVertLevels
+ theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
+ tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ end do
- end if
end do
end if
@@ -856,14 +846,12 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
- end do
+ do k=1,grid % nVertLevels
+ delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ end do
- end if
end do
do iCell = 1, nCells
@@ -876,17 +864,15 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- theta_turb_flux = h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * theta_turb_flux
+ do k=1,grid % nVertLevels
+ theta_turb_flux = h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
- end do
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
- end if
end do
deallocate(delsq_theta)
@@ -903,14 +889,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) )
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
- end do
- end if
+ do k=1,grid % nVertLevels
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) )
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
end do
else if (config_theta_adv_order == 3) then
@@ -918,37 +902,33 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
+ do k=1,grid % nVertLevels
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(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) * theta(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) * theta(k,grid % CellsOnCell % array (i,cell2))
- end do
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ end do
! 3rd order stencil
- if( u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
- else
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
- end if
-
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
-
- end do
- end if
+ if( u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+ else
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+ end if
+
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+
+ end do
end do
else if (config_theta_adv_order == 4) then
@@ -956,30 +936,25 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
+ do k=1,grid % nVertLevels
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(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) * theta(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) * theta(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
- end do
-
- end if
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
end do
end if
@@ -1162,18 +1137,16 @@
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
- u(k,iEdge) = u(k,iEdge) + dt*tend_u(k,iEdge) &
- -(0.5*dt/dcEdge(iEdge))*( &
- (geopotential(k+1,cell2)-geopotential(k+1,cell1)) &
- +(geopotential(k ,cell2)-geopotential(k ,cell1)) &
- +cqu(k,iEdge)*(alpha(k,cell2)+alpha(k,cell1))* &
- 0.5*(pressure(k+1,cell2)-pressure(k+1,cell1) &
- +pressure(k ,cell2)-pressure(k ,cell1))) &
- -smext*dcEdge(iEdge)*(dpsdt(cell2)-dpsdt(cell1))/h_edge(k,iEdge)
- end do
- end if
+ do k=1,nVertLevels
+ u(k,iEdge) = u(k,iEdge) + dt*tend_u(k,iEdge) &
+ -(0.5*dt/dcEdge(iEdge))*( &
+ (geopotential(k+1,cell2)-geopotential(k+1,cell1)) &
+ +(geopotential(k ,cell2)-geopotential(k ,cell1)) &
+ +cqu(k,iEdge)*(alpha(k,cell2)+alpha(k,cell1))* &
+ 0.5*(pressure(k+1,cell2)-pressure(k+1,cell1) &
+ +pressure(k ,cell2)-pressure(k ,cell1))) &
+ -smext*dcEdge(iEdge)*(dpsdt(cell2)-dpsdt(cell1))/h_edge(k,iEdge)
+ end do
end do
!
@@ -1184,14 +1157,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) + flux
- tend_h(k,cell2) = tend_h(k,cell2) - flux
- end do
- end if
do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) + flux
+ tend_h(k,cell2) = tend_h(k,cell2) - flux
+ end do
+ do k=1,nVertLevels
uhAvg(k,iEdge) = uhAvg(k,iEdge) + u(k,iEdge) * h_edge(k,iEdge)
end do
end do
@@ -1237,18 +1208,16 @@
!
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5*(h(k,cell1)+h(k,cell2)) ! here is update of h_edge
- he_old = 0.5*(h_old(k,cell1)+h_old(k,cell2))
- flux = 0.5*(u(k,iEdge) * h_edge(k,iEdge) - u_old(k,iEdge) * he_old)* &
- (theta_old(k,cell1)+theta_old(k,cell2))*dvEdge(iEdge)
- theta(k,cell1) = theta(k,cell1) - dt*flux/areaCell(cell1)
- theta(k,cell2) = theta(k,cell2) + dt*flux/areaCell(cell2)
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = 0.5*(h(k,cell1)+h(k,cell2)) ! here is update of h_edge
+ he_old = 0.5*(h_old(k,cell1)+h_old(k,cell2))
+ flux = 0.5*(u(k,iEdge) * h_edge(k,iEdge) - u_old(k,iEdge) * he_old)* &
+ (theta_old(k,cell1)+theta_old(k,cell2))*dvEdge(iEdge)
+ theta(k,cell1) = theta(k,cell1) - dt*flux/areaCell(cell1)
+ theta(k,cell2) = theta(k,cell2) + dt*flux/areaCell(cell2)
+ end do
end do
@@ -1360,16 +1329,14 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- 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
+ do k=1,grid % nVertLevels
+ 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 if
+ end do
end do
else if (config_scalar_adv_order == 3) then
@@ -1377,53 +1344,49 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
+ do k=1,grid % nVertLevels
- 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(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(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(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(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
+ 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)
+ d2fdx2_cell1 = d2fdx2_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)
+ d2fdx2_cell2 = d2fdx2_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(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(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
+
! 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(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(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
-! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
-! end if
-
- 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
+! if (uhAvg(k,iEdge) > 0) then
+! flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+! 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(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+! end if
+
+ 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 if
+ end do
end do
else if (config_scalar_adv_order == 4) then
@@ -1431,33 +1394,29 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
+ do k=1,grid % nVertLevels
- 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(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(iScalar,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
- 0.5*(scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- 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
+ 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)
+ d2fdx2_cell1 = d2fdx2_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)
+ d2fdx2_cell2 = d2fdx2_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(iScalar,k,cell1) + scalar_new(iScalar,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ 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 if
+ end do
end do
end if
@@ -1604,19 +1563,17 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- cell_upwind = cell2
- if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
- 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
+ cell_upwind = cell2
+ if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
+ 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 do
else if (config_scalar_adv_order >= 3) then
@@ -1624,44 +1581,40 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- cell_upwind = cell2
- if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
- do iScalar=1,num_scalars
+ cell_upwind = cell2
+ if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
+ 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)
+ d2fdx2_cell1 = d2fdx2_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)
+ d2fdx2_cell2 = d2fdx2_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(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(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
- 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(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(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(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(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(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
+ 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 do
end if
@@ -1691,24 +1644,22 @@
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 iScalar=1,num_scalars
+ do iScalar=1,num_scalars
- 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))
+ 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
- fdir = 1.0
- else
- fdir = -1.0
- end if
- 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
+ iEdge = grid % EdgesOnCell % array (i,iCell)
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ fdir = 1.0
+ else
+ fdir = -1.0
+ end if
+ 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 do
@@ -1747,17 +1698,15 @@
do iEdge = 1, grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do iScalar=1,num_scalars
- flux = h_flux(iScalar,iEdge)
- if (flux > 0) then
- flux = flux * min(scale_out(iScalar,cell1,km0), scale_in(iScalar,cell2,km0))
- else
- flux = flux * min(scale_in(iScalar,cell1,km0), scale_out(iScalar,cell2,km0))
- end if
- h_flux(iScalar,iEdge) = flux
- end do
- end if
+ do iScalar=1,num_scalars
+ flux = h_flux(iScalar,iEdge)
+ if (flux > 0) then
+ flux = flux * min(scale_out(iScalar,cell1,km0), scale_in(iScalar,cell2,km0))
+ else
+ flux = flux * min(scale_in(iScalar,cell1,km0), scale_out(iScalar,cell2,km0))
+ end if
+ h_flux(iScalar,iEdge) = flux
+ end do
end do
! rescale the vertical flux
@@ -1792,14 +1741,12 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- 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
+ 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 do
! decouple from mass
@@ -1906,11 +1853,9 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
end do
!
@@ -1918,16 +1863,10 @@
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
- if (verticesOnEdge(2,iEdge) > 0) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
+ do k=1,nVertLevels
+ circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
+ circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
+ end do
end do
do iVertex=1,nVertices
do k=1,nVertLevels
@@ -1943,16 +1882,10 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- end do
- end if
- if(cell2 > 0) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- end do
- end if
+ do k=1,nVertLevels
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ end do
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
@@ -1985,11 +1918,9 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
- do k = 1,nVertLevels
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end if
+ do k = 1,nVertLevels
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
end do
end do
@@ -2001,7 +1932,7 @@
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -2036,12 +1967,10 @@
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
- do k=1,nVertLevels
+ iEdge = edgesOnVertex(i,iVertex)
+ do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- end do
- end if
+ end do
end do
end do
! tdr
@@ -2065,12 +1994,10 @@
pv_cell(:,:) = 0.0
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
- do k = 1,nVertLevels
+ 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)
- end do
- end if
+ end do
end do
end do
! tdr
@@ -2082,12 +2009,10 @@
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
- do k = 1,nVertLevels
+ do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
- end do
- end if
+ end do
end do
! tdr
@@ -2102,4 +2027,84 @@
end subroutine compute_solve_diagnostics
+
+ subroutine compute_w (s_new, s_old, grid, dt )
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute (diagnose) vertical velocity (used by physics)
+ !
+ ! Input: s_new - current model state
+ ! s_old - previous model state
+ ! grid - grid metadata
+ ! dt - timestep between new and old
+ !
+ ! Output: w (m/s)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: s_new
+ type (grid_state), intent(in) :: s_old
+ type (grid_meta), intent(inout) :: grid
+
+ real (kind=RKIND), intent(in) :: dt
+
+ real (kind=RKIND), dimension(:,:), pointer :: geo_new, geo_old, u, ww, h, w
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, rdnw, fnm, fnp
+
+ integer :: iEdge, iCell, k, cell1, cell2
+ real (kind=RKIND), dimension( grid % nVertlevels + 1 ) :: wdwn
+ real (kind=RKIND) :: flux
+
+ geo_new => s_new % geopotential % array
+ geo_old => s_old % geopotential % array
+ u => s_new % u % array
+ ww => s_new % ww % array
+ h => s_new % h % array
+ w => s_new % w % array
+ dvEdge => grid % dvEdge % array
+ rdnw => grid % rdnw % array
+ fnm => grid % fnm % array
+ fnp => grid % fnp % array
+
+ w = 0.
+
+ do iCell=1, grid % nCellsSolve
+ wdwn(1) = 0.
+ do k=2,grid % nVertlevels + 1
+ wdwn(k) = (0.5*(ww(k,iCell)+ww(k-1,iCell))/h(k-1,iCell)) &
+ *rdnw(k-1)*(geo_new(k,iCell)-geo_new(k-1,iCell))
+ enddo
+ w(1,iCell) = 0.
+ do k=2, grid % nVertLevels
+ w(k,iCell) = (((geo_new(k,iCell)-geo_old(k,iCell))/dt)+ &
+ fnm(k)*wdwn(k+1)+fnp(k)*wdwn(k))/gravity
+ enddo
+ k = grid % nVertLevels + 1
+ w(k,iCell) = ((geo_new(k,iCell)-geo_old(k,iCell))/dt)/gravity
+ enddo
+
+ 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 k=2, grid % nVertLevels
+ flux = 0.25*(u(k,iEdge)+u(k-1,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+ enddo
+ k = 1
+ flux = 0.5*(1.5*u(1,iEdge)-0.5*u(2,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+
+ k = grid % nVertLevels + 1
+ flux = 0.5*(1.5*u(k-1,iEdge)-0.5*u(k-2,iEdge)*(geo_new(k,cell2)-geo_new(k,cell1))) * dvEdge(iEdge) / gravity
+ w(k,cell1) = w(k,cell1) + flux/ grid % AreaCell % array (cell1)
+ w(k,cell2) = w(k,cell2) + flux/ grid % AreaCell % array (cell2)
+
+ end if
+ end do
+
+ end subroutine compute_w
+
end module time_integration
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_grid_types.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_grid_types.F 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_grid_types.F 2010-06-01 20:30:27 UTC (rev 320)
@@ -70,6 +70,9 @@
#include "field_dimensions.inc"
+ logical :: on_a_sphere
+ real (kind=RKIND) :: sphere_radius
+
#include "time_invariant_fields.inc"
end type grid_meta
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_input.F 2010-06-01 20:30:27 UTC (rev 320)
@@ -49,6 +49,9 @@
integer :: i, j, k
type (io_input_object) :: input_obj
#include "dim_decls.inc"
+
+ character (len=16) :: c_on_a_sphere
+ real (kind=RKIND) :: r_sphere_radius
integer :: readCellStart, readCellEnd, nReadCells
integer :: readEdgeStart, readEdgeEnd, nReadEdges
@@ -729,6 +732,18 @@
#include "dim_dummy_args.inc"
)
+ !
+ ! Read attributes
+ !
+ call io_input_get_att_text(input_obj, 'on_a_sphere', c_on_a_sphere)
+ call io_input_get_att_real(input_obj, 'sphere_radius', r_sphere_radius)
+ if (index(c_on_a_sphere, 'YES') /= 0) then
+ domain % blocklist % mesh % on_a_sphere = .true.
+ else
+ domain % blocklist % mesh % on_a_sphere = .false.
+ end if
+ domain % blocklist % mesh % sphere_radius = r_sphere_radius
+
if (.not. config_do_restart) then
input_obj % time = 1
else
@@ -1100,7 +1115,59 @@
end subroutine io_input_get_dimension
+
+ subroutine io_input_get_att_real(input_obj, attname, attvalue)
+
+ implicit none
+ type (io_input_object), intent(in) :: input_obj
+ character (len=*), intent(in) :: attname
+ real (kind=RKIND), intent(out) :: attvalue
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ if (RKIND == 8) then
+ nferr = nf_get_att_double(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
+ else
+ nferr = nf_get_att_real(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
+ end if
+ if (nferr /= NF_NOERR) then
+ write(0,*) 'Warning: Attribute '//trim(attname)//' not found in '//trim(input_obj % filename)
+ if (index(attname, 'sphere_radius') /= 0) then
+ write(0,*) ' Setting '//trim(attname)//' to 1.0'
+ attvalue = 1.0
+ end if
+ end if
+
+ end subroutine io_input_get_att_real
+
+
+ subroutine io_input_get_att_text(input_obj, attname, attvalue)
+
+ implicit none
+
+ type (io_input_object), intent(in) :: input_obj
+ character (len=*), intent(in) :: attname
+ character (len=*), intent(out) :: attvalue
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ nferr = nf_get_att_text(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
+ if (nferr /= NF_NOERR) then
+ write(0,*) 'Warning: Attribute '//trim(attname)//' not found in '//trim(input_obj % filename)
+ if (index(attname, 'on_a_sphere') /= 0) then
+ write(0,*) ' Setting '//trim(attname)//' to ''YES'''
+ attvalue = 'YES'
+ end if
+ end if
+
+ end subroutine io_input_get_att_text
+
+
subroutine io_input_field0dReal(input_obj, field)
implicit none
Modified: branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/framework/module_io_output.F 2010-06-01 20:30:27 UTC (rev 320)
@@ -90,6 +90,7 @@
! although in future, work needs to be done to write model state
! from many distributed blocks
call io_output_init(output_obj, domain % dminfo, &
+ block_ptr % mesh, &
#include "output_dim_actual_args.inc"
)
@@ -326,6 +327,7 @@
subroutine io_output_init( output_obj, &
dminfo, &
+ mesh, &
#include "dim_dummy_args.inc"
)
@@ -335,6 +337,7 @@
type (io_output_object), intent(inout) :: output_obj
type (dm_info), intent(in) :: dminfo
+ type (grid_meta), intent(in) :: mesh
#include "dim_dummy_decls.inc"
integer :: nferr
@@ -348,6 +351,17 @@
#endif
#include "netcdf_def_dims_vars.inc"
+
+ if (mesh % on_a_sphere) then
+ nferr = nf_put_att_text(output_obj % wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, 'YES ')
+ else
+ nferr = nf_put_att_text(output_obj % wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, 'NO ')
+ end if
+ if (RKIND == 8) then
+ nferr = nf_put_att_double(output_obj % wr_ncid, NF_GLOBAL, 'sphere_radius', NF_DOUBLE, 1, mesh % sphere_radius)
+ else
+ nferr = nf_put_att_real(output_obj % wr_ncid, NF_GLOBAL, 'sphere_radius', NF_FLOAT, 1, mesh % sphere_radius)
+ end if
nferr = nf_enddef(output_obj % wr_ncid)
end if
Modified: branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/registry/gen_inc.c 2010-06-01 20:30:27 UTC (rev 320)
@@ -222,12 +222,14 @@
fd = fopen("field_dimensions.inc", "w");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="red">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sSolve</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sSolve</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sSolve</font>
<font color="gray">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
@@ -239,12 +241,17 @@
*/
fd = fopen("dim_dummy_args.inc", "w");
dim_ptr = dims;
- if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
+ if (dim_ptr && dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) {
fortprintf(fd, " %s", dim_ptr->name_in_code);
dim_ptr = dim_ptr->next;
}
+ else if (dim_ptr && dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) {
+ fortprintf(fd, " %s", dim_ptr->name_in_file);
+ dim_ptr = dim_ptr->next;
+ }
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, " &</font>
<font color="gray">");
@@ -257,12 +264,17 @@
*/
fd = fopen("dim_dummy_decls.inc", "w");
dim_ptr = dims;
- if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
+ if (dim_ptr && dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) {
fortprintf(fd, " integer, intent(in) :: %s", dim_ptr->name_in_code);
dim_ptr = dim_ptr->next;
}
+ else if (dim_ptr && dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) {
+ fortprintf(fd, " integer, intent(in) :: %s", dim_ptr->name_in_file);
+ dim_ptr = dim_ptr->next;
+ }
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -276,7 +288,8 @@
fd = fopen("dim_decls.inc", "w");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %s</font>
<font color="gray">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
@@ -289,7 +302,8 @@
fd = fopen("read_dims.inc", "w");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " call io_input_get_dimension(input_obj, \'%s\', %s)</font>
<font color="blue">", dim_ptr->name_in_file, dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " call io_input_get_dimension(input_obj, \'%s\', %s)</font>
<font color="blue">", dim_ptr->name_in_file, dim_ptr->name_in_code);
+ else if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " call io_input_get_dimension(input_obj, \'%s\', %s)</font>
<font color="gray">", dim_ptr->name_in_file, dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
@@ -365,7 +379,8 @@
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " g %% %s = %s</font>
<font color="blue">", dim_ptr->name_in_code, dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " g %% %s = %s</font>
<font color="blue">", dim_ptr->name_in_code, dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " g %% %s = %s</font>
<font color="black">", dim_ptr->name_in_file, dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -392,7 +407,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
@@ -400,7 +416,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -432,7 +449,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
@@ -440,7 +458,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -541,7 +560,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
@@ -549,7 +569,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -582,7 +603,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
}
else {
fortprintf(fd, "%i", dimlist_ptr->dim->constant_value);
@@ -595,7 +617,8 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
else
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
}
else {
fortprintf(fd, ", %i", dimlist_ptr->dim->constant_value);
@@ -783,7 +806,7 @@
fortprintf(fd, " integer :: rdDimIDTime</font>
<font color="red">");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: rdDimID%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: rdDimID%s</font>
<font color="black">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -791,7 +814,7 @@
fortprintf(fd, " integer :: rdLocalTime</font>
<font color="red">");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: rdLocal%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: rdLocal%s</font>
<font color="black">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -851,7 +874,8 @@
fortprintf(fd, "#endif</font>
<font color="red">");
}
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " %s%id %% ioinfo %% count(%i) = block %% mesh %% %s</font>
<font color="blue">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " %s%id %% ioinfo %% count(%i) = block %% mesh %% %s</font>
<font color="blue">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " %s%id %% ioinfo %% count(%i) = block %% mesh %% %s</font>
<font color="black">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " %s%id %% ioinfo %% count(%i) = %s</font>
<font color="gray">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
if (has_vert_dim) {
@@ -884,7 +908,8 @@
if (i < var_ptr->ndims) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
}
@@ -904,7 +929,8 @@
while (dimlist_ptr) {
if (i < var_ptr->ndims) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
}
@@ -930,7 +956,8 @@
if (i < var_ptr->ndims) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
}
@@ -938,7 +965,8 @@
i++;
while (dimlist_ptr) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
@@ -994,7 +1022,8 @@
if (i < var_ptr->ndims)
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " %s", dimlist_ptr->dim->name_in_code);
else {
@@ -1014,7 +1043,8 @@
while (dimlist_ptr) {
if (i < var_ptr->ndims)
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
else {
@@ -1046,7 +1076,8 @@
dimlist_ptr = var_ptr->dimlist;
while (i <= var_ptr->ndims) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " do i%i=1,block %% mesh %% %s</font>
<font color="blue">", i, dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " do i%i=1,block %% mesh %% %s</font>
<font color="blue">", i, dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " do i%i=1,block %% mesh %% %s</font>
<font color="black">", i, dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " do i%i=1,%s</font>
<font color="gray">", i, dimlist_ptr->dim->name_in_code);
@@ -1104,7 +1135,7 @@
fortprintf(fd, " nferr = nf_inq_dimlen(input_obj %% rd_ncid, input_obj %% rdDimIDTime, input_obj %% rdLocalTime)</font>
<font color="red">");
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) {
fortprintf(fd, " nferr = nf_inq_dimid(input_obj %% rd_ncid, \'%s\', input_obj %% rdDimID%s)</font>
<font color="black">", dim_ptr->name_in_file, dim_ptr->name_in_file);
fortprintf(fd, " nferr = nf_inq_dimlen(input_obj %% rd_ncid, input_obj %% rdDimID%s, input_obj %% rdLocal%s)</font>
<font color="gray">", dim_ptr->name_in_file, dim_ptr->name_in_file);
}
@@ -1128,13 +1159,25 @@
dim_ptr = dims;
while (dim_ptr->constant_value >= 0 || is_derived_dim(dim_ptr->name_in_code)) dim_ptr = dim_ptr->next;
- fortprintf(fd, " if (trim(dimname) == \'%s\') then</font>
<font color="red">", dim_ptr->name_in_code);
- fortprintf(fd, " dimsize = input_obj %% rdLocal%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ if (!dim_ptr->namelist_defined) {
+ fortprintf(fd, " if (trim(dimname) == \'%s\') then</font>
<font color="blue">", dim_ptr->name_in_code);
+ fortprintf(fd, " dimsize = input_obj %% rdLocal%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ }
+ else {
+ fortprintf(fd, " if (trim(dimname) == \'%s\') then</font>
<font color="blue">", dim_ptr->name_in_file);
+ fortprintf(fd, " dimsize = %s</font>
<font color="red">", dim_ptr->name_in_code);
+ }
dim_ptr = dim_ptr->next;
while (dim_ptr) {
if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
- fortprintf(fd, " else if (trim(dimname) == \'%s\') then</font>
<font color="red">", dim_ptr->name_in_code);
- fortprintf(fd, " dimsize = input_obj %% rdLocal%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ if (!dim_ptr->namelist_defined) {
+ fortprintf(fd, " else if (trim(dimname) == \'%s\') then</font>
<font color="blue">", dim_ptr->name_in_code);
+ fortprintf(fd, " dimsize = input_obj %% rdLocal%s</font>
<font color="blue">", dim_ptr->name_in_file);
+ }
+ else {
+ fortprintf(fd, " else if (trim(dimname) == \'%s\') then</font>
<font color="blue">", dim_ptr->name_in_file);
+ fortprintf(fd, " dimsize = %s</font>
<font color="gray">", dim_ptr->name_in_code);
+ }
}
dim_ptr = dim_ptr->next;
}
@@ -1252,7 +1295,8 @@
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sGlobal</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sGlobal</font>
<font color="blue">", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " integer :: %sGlobal</font>
<font color="black">", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -1267,7 +1311,8 @@
dim_ptr = dims;
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " %sGlobal = block_ptr %% mesh %% %s</font>
<font color="blue">", dim_ptr->name_in_code, dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " %sGlobal = block_ptr %% mesh %% %s</font>
<font color="blue">", dim_ptr->name_in_code, dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, " %sGlobal = block_ptr %% mesh %% %s</font>
<font color="black">", dim_ptr->name_in_file, dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, "</font>
<font color="gray">");
@@ -1281,11 +1326,13 @@
fd = fopen("output_dim_actual_args.inc", "w");
dim_ptr = dims;
if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
- fortprintf(fd, " %sGlobal", dim_ptr->name_in_code);
+ if (!dim_ptr->namelist_defined) fortprintf(fd, " %sGlobal", dim_ptr->name_in_code);
+ else fortprintf(fd, " %sGlobal", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
while (dim_ptr) {
- if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %sGlobal", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && !dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %sGlobal", dim_ptr->name_in_code);
+ if (dim_ptr->constant_value < 0 && dim_ptr->namelist_defined && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %sGlobal", dim_ptr->name_in_file);
dim_ptr = dim_ptr->next;
}
fortprintf(fd, " &</font>
<font color="gray">");
@@ -1393,7 +1440,8 @@
if (i < var_ptr->ndims) {
fortprintf(fd, " %s%id %% ioinfo %% start(%i) = 1</font>
<font color="red">", vtype, var_ptr->ndims, i);
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " %s%id %% ioinfo %% count(%i) = domain %% blocklist %% mesh %% %s</font>
<font color="blue">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " %s%id %% ioinfo %% count(%i) = domain %% blocklist %% mesh %% %s</font>
<font color="blue">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " %s%id %% ioinfo %% count(%i) = domain %% blocklist %% mesh %% %s</font>
<font color="black">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " %s%id %% ioinfo %% count(%i) = %s</font>
<font color="gray">", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
}
@@ -1418,7 +1466,8 @@
if (i < var_ptr->ndims)
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
else {
@@ -1437,7 +1486,8 @@
while (dimlist_ptr) {
if (i < var_ptr->ndims)
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
else {
@@ -1463,7 +1513,8 @@
dimlist_ptr = var_ptr->dimlist;
while (dimlist_ptr) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
@@ -1480,7 +1531,8 @@
dimlist_ptr = var_ptr->dimlist;
while (i <= var_ptr->ndims) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " do i%i=1,domain %% blocklist %% mesh %% %s</font>
<font color="blue">", i, dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " do i%i=1,domain %% blocklist %% mesh %% %s</font>
<font color="blue">", i, dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " do i%i=1,domain %% blocklist %% mesh %% %s</font>
<font color="black">", i, dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " do i%i=1,%s</font>
<font color="gray">", i, dimlist_ptr->dim->name_in_code);
@@ -1525,7 +1577,8 @@
dimlist_ptr = var_ptr->dimlist;
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, " domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, " domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, " domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, " %s", dimlist_ptr->dim->name_in_code);
@@ -1533,7 +1586,8 @@
i++;
while (dimlist_ptr) {
if (dimlist_ptr->dim->constant_value < 0)
- fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_file);
else
fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
Modified: branches/ocean_projects/z_level_mrp/mpas/src/registry/parse.c
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/registry/parse.c 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/registry/parse.c 2010-06-01 20:30:27 UTC (rev 320)
@@ -49,6 +49,7 @@
{
char word[1024];
struct namelist * nls_ptr;
+ struct namelist * nls_chk_ptr;
struct dimension * dim_ptr;
struct variable * var_ptr;
struct dimension_list * dimlist_ptr;
@@ -96,9 +97,31 @@
else if (strncmp(word, "dim", 1024) == 0) {
NEW_DIMENSION(dim_ptr->next)
dim_ptr = dim_ptr->next;
+ dim_ptr->namelist_defined = 0;
getword(regfile, dim_ptr->name_in_file);
getword(regfile, dim_ptr->name_in_code);
dim_ptr->constant_value = is_integer_constant(dim_ptr->name_in_code);
+ if (strncmp(dim_ptr->name_in_code, "namelist:", 9) == 0) {
+ dim_ptr->namelist_defined = 1;
+ sprintf(dim_ptr->name_in_code, "%s", (dim_ptr->name_in_code)+9);
+
+ /* Check that the referenced namelist variable is defined as an integer variable */
+ nls_chk_ptr = (*nls)->next;
+ while (nls_chk_ptr) {
+ if (strncmp(nls_chk_ptr->name, dim_ptr->name_in_code, 1024) == 0) {
+ if (nls_chk_ptr->vtype != INTEGER) {
+ printf("</font>
<font color="black">Registry error: Namelist variable %s must be an integer for namelist-derived dimension %s</font>
<font color="black"></font>
<font color="blue">", nls_chk_ptr->name, dim_ptr->name_in_file);
+ return 1;
+ }
+ break;
+ }
+ nls_chk_ptr = nls_chk_ptr->next;
+ }
+ if (!nls_chk_ptr) {
+ printf("</font>
<font color="black">Registry error: Namelist variable %s not defined for namelist-derived dimension %s</font>
<font color="black"></font>
<font color="gray">", dim_ptr->name_in_code, dim_ptr->name_in_file);
+ return 1;
+ }
+ }
}
else if (strncmp(word, "var", 1024) == 0) {
NEW_VARIABLE(var_ptr->next)
Modified: branches/ocean_projects/z_level_mrp/mpas/src/registry/registry_types.h
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/registry/registry_types.h 2010-06-01 17:17:03 UTC (rev 319)
+++ branches/ocean_projects/z_level_mrp/mpas/src/registry/registry_types.h 2010-06-01 20:30:27 UTC (rev 320)
@@ -31,6 +31,7 @@
char name_in_file[1024];
char name_in_code[1024];
int constant_value;
+ int namelist_defined;
struct dimension * next;
};
</font>
</pre>