<p><b>laura@ucar.edu</b> 2010-05-20 15:18:40 -0600 (Thu, 20 May 2010)</p><p>Final updates between trunk and branch atmos_physics<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/namelist.input.ocean
===================================================================
--- branches/atmos_physics/namelist.input.ocean 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/namelist.input.ocean 2010-05-20 21:18:40 UTC (rev 295)
@@ -1,21 +1,21 @@
&sw_model
- config_test_case = 5
+ config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 300.0
- config_ntimesteps = 3000
- config_output_interval = 300
- config_stats_interval = 10
- config_visc = 4.0
+ config_dt = 120.0
+ config_ntimesteps = 30
+ config_output_interval = 30
+ config_stats_interval = 10
+ config_visc = 100.0
/
&io
- config_input_name = grid.nc
- config_output_name = output.nc
- config_restart_name = restart.nc
+ config_input_name = 'grid.quad.20km.nc'
+ config_output_name = 'output.quad.20km.nc'
+ config_restart_name = 'restart.nc'
/
&restart
- config_restart_interval = 864000
+ config_restart_interval = 3000
config_do_restart = .false.
config_restart_time = 1036800.0
/
Modified: branches/atmos_physics/src/core_hyd_atmos/Makefile
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -19,7 +19,7 @@
mpas_interface.o: module_advection.o module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
Modified: branches/atmos_physics/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_advection.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_hyd_atmos/module_advection.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -96,7 +96,7 @@
do_the_cell = .true.
do i=1,n
- if (cell_list(i) <= 0) do_the_cell = .false.
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
end do
Modified: branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -466,9 +466,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
@@ -755,33 +753,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)
@@ -794,16 +784,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)
@@ -904,16 +888,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
@@ -927,14 +909,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
@@ -947,17 +927,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)
@@ -974,14 +952,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
@@ -989,37 +965,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
@@ -1027,30 +999,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
@@ -1233,18 +1200,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
!
@@ -1255,14 +1220,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
@@ -1308,18 +1271,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
@@ -1431,16 +1392,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
@@ -1448,53 +1407,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
@@ -1502,33 +1457,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
@@ -1675,19 +1626,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
@@ -1695,44 +1644,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
-
- 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
+ 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
- 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
+ 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 do
end if
@@ -1762,24 +1707,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
@@ -1818,17 +1761,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
@@ -1863,14 +1804,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
@@ -1977,11 +1916,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
!
@@ -1989,16 +1926,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
@@ -2014,16 +1945,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)
@@ -2056,11 +1981,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
@@ -2072,7 +1995,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
@@ -2107,12 +2030,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
@@ -2136,12 +2057,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
@@ -2153,12 +2072,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
Modified: branches/atmos_physics/src/core_ocean/Makefile
===================================================================
--- branches/atmos_physics/src/core_ocean/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_ocean/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -19,10 +19,9 @@
mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
Modified: branches/atmos_physics/src/core_ocean/Registry
===================================================================
--- branches/atmos_physics/src/core_ocean/Registry 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_ocean/Registry 2010-05-20 21:18:40 UTC (rev 295)
@@ -6,9 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
-# mrp 100120:
namelist integer sw_model config_stats_interval 100
-# mrp 100120 end
namelist real sw_model config_visc 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -85,7 +83,8 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
var real u_src ( nVertLevels nEdges ) iro u_src - -
# Prognostic variables: read from input, saved in restart, and written to output
@@ -105,14 +104,15 @@
var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
-# mrp 100112:
-var real zmid ( nVertLevels nCells Time ) o zmid - -
-var real zbot ( nVertLevels nCells Time ) o zbot - -
+var real zMid ( nVertLevels nCells Time ) o zMid - -
+var real zBot ( nVertLevels nCells Time ) o zBot - -
var real zSurface ( nCells Time ) o zSurface - -
+var real zMidEdge ( nVertLevels nEdges Time ) o zMidEdge - -
+var real zBotEdge ( nVertLevels nEdges Time ) o zBotEdge - -
+var real zSurfaceEdge ( nEdges Time ) o zSurfaceEdge - -
var real pmid ( nVertLevels nCells Time ) o pmid - -
var real pbot ( nVertLevels nCells Time ) o pbot - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
-# mrp 100112 end
# Other diagnostic variables: neither read nor written to any files
Modified: branches/atmos_physics/src/core_ocean/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_ocean/module_test_cases.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_ocean/module_test_cases.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -23,9 +23,6 @@
integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND) :: delta_rho
- integer :: nCells, nEdges, nVertices, nVertLevels
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -76,84 +73,16 @@
stop
end if
- ! mrp 100121:
- !
- ! Initialize u_src, rho, alpha
- ! This is a temporary fix until everything is specified correctly
- ! in an initial conditions file.
- !
block_ptr => domain % blocklist
do while (associated(block_ptr))
- h => block_ptr % time_levs(1) % state % h % array
- u => block_ptr % time_levs(1) % state % u % array
- rho => block_ptr % time_levs(1) % state % rho % array
- u_src => block_ptr % mesh % u_src % array
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % time_levs(1) % state, &
+ block_ptr % time_levs(i) % state)
+ end do
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
- ! Momentum forcing u_src
- if (config_test_case > 0) then
- ! for shallow water test cases:
- u_src = 0.0
- elseif (config_test_case == 0 ) then
- ! for rectangular basin:
- do iEdge=1,nEdges
- u_src(1:nVertLevels, iEdge) = u_src(1, iEdge) / nVertLevels
- end do
- endif
-
- ! define the density for multiple layers
- delta_rho=0.0
- do iLevel = 1,nVertLevels
- rho(iLevel,1) = delta_rho*(iLevel-1)
- enddo
- rho(:,1) = rho(:,1) - sum(rho(:,1))/nVertLevels + 1000.0
- do iLevel = 1,nVertLevels
- rho(iLevel,:) = rho(iLevel,1)
- enddo
-
-#ifdef EXPAND_LEVELS
- print '(10a)', 'EXPAND_LEVELS compiler flag is on.', &
- ' Copying h and u from k=1 to other levels.'
- print '(a,i5)', 'EXPAND_LEVELS =', EXPAND_LEVELS
- print '(a,i5)', 'nVertLevels =', nVertLevels
- do iCell=1,nCells
- ! make the total thickness equal to the single-layer thickness:
- h(1:nVertLevels, iCell) = h(1,iCell) / nVertLevels
- end do
-
- do iEdge=1,nEdges
- u(2:nVertLevels, iEdge) = u(1,iEdge)
- end do
-#else
- print '(10a)', 'EXPAND_LEVELS compiler flag is off.'
-#endif
-
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, &
- block_ptr % time_levs(i) % state)
- end do
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min h max h ',&
- ' min u_src max u_src '
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,:)), maxval(u(iLevel,:)), &
- minval(h(iLevel,:)), maxval(h(iLevel,:)), &
- minval(u_src(iLevel,:)), maxval(u_src(iLevel,:))
- enddo
-
- block_ptr => block_ptr % next
+ block_ptr => block_ptr % next
end do
- ! mrp 100121 end
end subroutine setup_sw_test_case
Modified: branches/atmos_physics/src/core_ocean/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_ocean/module_time_integration.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_ocean/module_time_integration.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -38,19 +38,18 @@
stop
end if
- block => domain % blocklist
- do while (associated(block))
- block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + dt
- ! mrp 100310 I added this to avoid running with NaNs
- if (isNaN(sum(block % time_levs(2) % state % u % array))) then
- print *, 'Stopping: NaN detected'
- call MPI_abort(MPI_COMM_WORLD,errorcode,ierr)
- endif
- ! mrp 100310 end
+ block => domain % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + dt
+ ! ! mrp 100310 I added this to avoid running with NaNs
+ ! if (isNaN(sum(block % time_levs(2) % state % u % array))) then
+ ! print *, 'Stopping: NaN detected'
+ ! call MPI_abort(MPI_COMM_WORLD,errorcode,ierr)
+ ! endif
+ ! ! mrp 100310 end
+ block => block % next
+ end do
- block => block % next
- end do
-
end subroutine timestep
@@ -136,7 +135,7 @@
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -253,12 +252,13 @@
type (grid_meta), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias, dist
+ real (kind=RKIND), allocatable, dimension(:) :: fluxVert
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zSurfaceEdge
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence
+ circulation, vorticity, ke, pv_edge, divergence, zBotEdge, zMidEdge
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: u_diffusion, visc
@@ -267,10 +267,8 @@
real (kind=RKIND), dimension(:,:), pointer :: MontPot
!mrp 100112 end
-!ocean
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
-!ocean
visc = config_visc
@@ -284,9 +282,10 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
vh => s % vh % array
- !mrp 100112:
MontPot => s % MontPot % array
- !mrp 100112 end
+ zBotEdge => s % zBotEdge % array
+ zSurfaceEdge=> s % zSurfaceEdge % array
+ zMidEdge => s % zMidEdge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -314,9 +313,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
-!ocean
u_src => grid % u_src % array
-!ocean
!
@@ -326,13 +323,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
end do
end if
- if (cell2 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -372,28 +369,45 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
- ! mrp 100112, this is the original shallow water formulation with grad H:
- !tend_u(k,iEdge) = &
- ! q &
- ! + u_diffusion &
- ! - ( ke(k,cell2) - ke(k,cell1) &
- ! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ! ) / dcEdge(iEdge)
- ! mrp 100112, changed to grad Montgomery potential:
+
tend_u(k,iEdge) = &
q &
+ u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) &
+ MontPot(k,cell2) - MontPot(k,cell1) &
) / dcEdge(iEdge)
- ! mrp 100112 end
+ end do
+ end do
-!ocean
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_src(k,iEdge)/rho_ref/h_edge(k,iEdge)
-!ocean
+ do iEdge=1,grid % nEdgesSolve
+ ! surface forcing
+ tend_u(1,iEdge) = tend_u(1,iEdge) + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- end do
- end do
+ ! bottom drag
+ tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) - u(nVertLevels,iEdge)/(100.0*86400.0)
+
+ enddo
+
+! vertical mixing
+ allocate(fluxVert(0:nVertLevels))
+ do iEdge=1,grid % nEdgesSolve
+ fluxVert(0) = 0.0
+ fluxVert(nVertLevels) = 0.0
+ do k=nVertLevels-1,1,-1
+ fluxVert(k) = ( u(k,iEdge) - u(k+1,iEdge) ) / (zMidEdge(k,iEdge) - zMidEdge(k+1,iEdge))
+ enddo
+ fluxVert = 1.0e-4 * fluxVert
+ do k=1,nVertLevels
+ if(k.eq.1) then
+ dist = zSurfaceEdge(iEdge) - zBotEdge(k,iEdge)
+ else
+ dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
+ endif
+ tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
+ enddo
+ enddo
+ deallocate(fluxVert)
+
#endif
#ifdef NCAR_FORMULATION
@@ -447,7 +461,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
@@ -493,15 +507,13 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- !mrp 100112:
real (kind=RKIND), dimension(:,:), pointer :: &
- zmid, zbot, pmid, pbot, MontPot, rho
- real (kind=RKIND), dimension(:), pointer :: zSurface
+ zMid, zBot, pmid, pbot, MontPot, rho, zBotEdge, zMidEdge
+ real (kind=RKIND), dimension(:), pointer :: zSurface, zSurfaceEdge
real (kind=RKIND) :: delta_p
character :: c1*6
- !mrp 100112 end
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -524,13 +536,15 @@
gradPVt => s % gradPVt % array
!mrp 100112:
rho => s % rho % array
- zmid => s % zmid % array
- zbot => s % zbot % array
+ zMid => s % zMid % array
+ zBot => s % zBot % array
+ zMidEdge => s % zMidEdge % array
+ zBotEdge => s % zBotEdge % array
+ zSurfaceEdge=> s % zSurfaceEdge % array
pmid => s % pmid % array
pbot => s % pbot % array
MontPot => s % MontPot % array
zSurface => s % zSurface % array
- !mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -555,7 +569,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -563,24 +577,39 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
+ elseif(cell1 <= nCells) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = h(k,cell1)
+ end do
+ elseif(cell2 <= nCells) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = h(k,cell2)
+ end do
end if
end do
+
!
+ ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
+ ! used to when reading for edges that do not exist
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) 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
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -592,6 +621,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -599,12 +629,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -640,7 +670,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -664,14 +694,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -679,14 +708,12 @@
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
-
+
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -698,7 +725,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -707,16 +733,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -727,7 +751,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -736,31 +759,29 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
+
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
@@ -770,25 +791,6 @@
enddo
!
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
-
- !mrp 100112:
- !
! Compute mid- and bottom-depth of each layer, from bottom up.
! Then compute mid- and bottom-pressure of each layer, and
! Montgomery Potential, from top down
@@ -798,14 +800,14 @@
! h_s is ocean topography: elevation above lowest point,
! and is positive. z coordinates are pos upward, and z=0
! is at lowest ocean point.
- zbot(nVertLevels,iCell) = h_s(iCell)
- zmid(nVertLevels,iCell) = zbot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
+ zBot(nVertLevels,iCell) = h_s(iCell)
+ zMid(nVertLevels,iCell) = zBot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
do k=nVertLevels-1,1,-1
- zbot(k,iCell) = zbot(k+1,iCell) + h(k+1,iCell)
- zmid(k,iCell) = zbot(k,iCell) + 0.5*h(k,iCell)
+ zBot(k,iCell) = zBot(k+1,iCell) + h(k+1,iCell)
+ zMid(k,iCell) = zBot(k,iCell) + 0.5*h(k,iCell)
end do
- ! rather than using zbot(0,iCell), I am adding another 2D variable.
- zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
+ ! rather than using zBot(0,iCell), I am adding another 2D variable.
+ zSurface(iCell) = zBot(1,iCell) + h(1,iCell)
! assume pressure at the surface is zero for now.
pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
@@ -821,18 +823,30 @@
+ pbot(k-1,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
end do
- !mrp 100112 end
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if(cell1<=nCells .and. cell2<=nCells) then
+ do k=1,nVertLevels
+ zBotEdge(k,iEdge) = (zBot(k,cell1)+zBot(k,cell2))/2.0
+ zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
+ enddo
+ zSurfaceEdge(iEdge) = (zSurface(cell1) + zSurface(cell2))/2.0
+ endif
+ enddo
+
+
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -841,7 +855,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -851,21 +865,21 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: branches/atmos_physics/src/core_sw/Makefile
===================================================================
--- branches/atmos_physics/src/core_sw/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_sw/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -16,10 +16,9 @@
mpas_interface.o: module_test_cases.o module_time_integration.o
clean:
- $(RM) *.o *.mod libdycore.a
+ $(RM) *.o *.mod *.f90 libdycore.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
- $(RM) $*.f90
Modified: branches/atmos_physics/src/core_sw/Registry
===================================================================
--- branches/atmos_physics/src/core_sw/Registry 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_sw/Registry 2010-05-20 21:18:40 UTC (rev 295)
@@ -81,7 +81,8 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
# Prognostic variables: read from input, saved in restart, and written to output
var real u ( nVertLevels nEdges Time ) iro u - -
Modified: branches/atmos_physics/src/core_sw/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_sw/module_time_integration.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/core_sw/module_time_integration.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -125,7 +125,7 @@
do while (associated(block))
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -297,13 +297,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
end do
end if
- if (cell2 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -343,6 +343,7 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
+
tend_u(k,iEdge) = &
q &
+ u_diffusion &
@@ -404,7 +405,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
@@ -450,7 +451,7 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -495,7 +496,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -503,7 +504,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
@@ -511,16 +512,22 @@
end do
!
+ ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
+ ! used to when reading for edges that do not exist
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) 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
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -532,6 +539,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -539,12 +547,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -580,7 +588,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -604,14 +612,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -623,10 +630,8 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -638,7 +643,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -647,16 +651,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -667,7 +669,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -676,30 +677,28 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
+
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
! Modify PV edge with upstream bias.
!
@@ -712,32 +711,38 @@
!
! set pv_edge = fEdge / h_edge at boundary points
!
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
+ ! if (maxval(boundaryEdge).ge.0) then
+ ! do iEdge = 1,nEdges
+ ! cell1 = cellsOnEdge(1,iEdge)
+ ! cell2 = cellsOnEdge(2,iEdge)
+ ! do k = 1,nVertLevels
+ ! if(boundaryEdge(k,iEdge).eq.1) then
+ ! v(k,iEdge) = 0.0
+ ! if(cell1.gt.0) then
+ ! h1 = h(k,cell1)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
+ ! h_edge(k,iEdge) = h1
+ ! else
+ ! h2 = h(k,cell2)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
+ ! h_edge(k,iEdge) = h2
+ ! endif
+ ! endif
+ ! enddo
+ ! enddo
+ ! endif
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -746,7 +751,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -756,22 +761,22 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
- tend_u => tend % u % array
+ boundaryEdge => grid % boundaryEdge % array
+ tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: branches/atmos_physics/src/driver/Makefile
===================================================================
--- branches/atmos_physics/src/driver/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/driver/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -10,10 +10,9 @@
mpas.o: module_subdriver.o
clean:
- $(RM) *.o *.mod
+ $(RM) *.o *.mod *.f90
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../core_$(CORE)
- $(RM) $*.f90
Modified: branches/atmos_physics/src/framework/Makefile
===================================================================
--- branches/atmos_physics/src/framework/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -35,13 +35,12 @@
module_io_output.o: module_grid_types.o module_dmpar.o module_sort.o module_configure.o
clean:
- $(RM) *.o *.mod libframework.a
+ $(RM) *.o *.mod *.f90 libframework.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES)
- $(RM) $*.f90
.c.o:
$(CC) $(CFLAGS) $(CPPFLAGS) $(CPPINCLUDES) -c $<
Modified: branches/atmos_physics/src/framework/module_block_decomp.F
===================================================================
--- branches/atmos_physics/src/framework/module_block_decomp.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/module_block_decomp.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -142,7 +142,11 @@
edgeIDListLocal(:) = edgeIDList(:)
do i=1,nEdges
- if (hash_search(h, cellsOnEdge(1,i))) then
+ do j=1,maxCells
+ if (cellsOnEdge(j,i) /= 0) exit
+ end do
+if (j > maxCells) write(0,*) 'Error in block_decomp_partitioned_edge_list: edge/vertex is not adjacent to any valid cells'
+ if (hash_search(h, cellsOnEdge(j,i))) then
lastEdge = lastEdge + 1
edgeIDList(lastEdge) = edgeIDListLocal(i)
else
@@ -231,8 +235,10 @@
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ end if
end if
end do
end do
@@ -264,10 +270,12 @@
write(0,*) 'block_decomp_add_halo: Somehow we don''t have the right number of non-ghost cells.'
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
- local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
- k = k + 1
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
+ k = k + 1
+ end if
end if
end do
end do
Modified: branches/atmos_physics/src/framework/module_dmpar.F
===================================================================
--- branches/atmos_physics/src/framework/module_dmpar.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/module_dmpar.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -707,8 +707,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- integer, dimension(nOwnedList), intent(in) :: arrayIn
- integer, dimension(nNeededList), intent(inout) :: arrayOut
+ integer, dimension(*), intent(in) :: arrayIn
+ integer, dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -784,7 +784,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -797,8 +797,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- integer, dimension(dim1,nOwnedList), intent(in) :: arrayIn
- integer, dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ integer, dimension(dim1,*), intent(in) :: arrayIn
+ integer, dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -876,7 +876,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -888,8 +888,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- real (kind=RKIND), dimension(nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -965,7 +965,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -978,8 +978,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1057,7 +1057,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -1070,8 +1070,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,dim2,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,dim2,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1149,7 +1149,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:,:) = arrayIn(:,:,:)
+ arrayOut(:,:,1:nNeededList) = arrayIn(:,:,1:nOwnedList)
end if
#endif
@@ -1161,7 +1161,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- integer, dimension(nField), intent(in) :: field
+ integer, dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1188,7 +1188,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- integer, dimension(ds:de,1:nField), intent(in) :: field
+ integer, dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1222,7 +1222,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(nField), intent(in) :: field
+ real (kind=RKIND), dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1249,7 +1249,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1283,7 +1283,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1321,7 +1321,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- integer, dimension(nField), intent(inout) :: field
+ integer, dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1348,7 +1348,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- integer, dimension(ds:de,1:nField), intent(inout) :: field
+ integer, dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1377,7 +1377,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(nField), intent(inout) :: field
+ real (kind=RKIND), dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1404,7 +1404,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1434,7 +1434,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1468,7 +1468,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1
- real (kind=RKIND), dimension(dim1), intent(inout) :: array
+ real (kind=RKIND), dimension(*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1528,7 +1528,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2
- real (kind=RKIND), dimension(dim1,dim2), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1592,7 +1592,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, dim3
- real (kind=RKIND), dimension(dim1,dim2,dim3), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
Modified: branches/atmos_physics/src/framework/module_io_input.F
===================================================================
--- branches/atmos_physics/src/framework/module_io_input.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/module_io_input.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -21,6 +21,7 @@
interface io_input_field
+ module procedure io_input_field0dReal
module procedure io_input_field1dReal
module procedure io_input_field2dReal
module procedure io_input_field3dReal
@@ -102,6 +103,7 @@
type (exchange_list), pointer :: sendEdgeList, recvEdgeList
type (exchange_list), pointer :: sendVertexList, recvVertexList
type (exchange_list), pointer :: sendVertLevelList, recvVertLevelList
+ type (exchange_list), pointer :: sendVertLevelP1List, recvVertLevelP1List
type (exchange_list), pointer :: send1Halo, recv1Halo
type (exchange_list), pointer :: send2Halo, recv2Halo
type (graph) :: partial_global_graph_info
@@ -692,7 +694,28 @@
deallocate(local_vertlevel_list)
deallocate(needed_vertlevel_list)
+ if (domain % dminfo % my_proc_id == 0) then
+ allocate(local_vertlevel_list(nVertLevels+1))
+ do i=1,nVertLevels+1
+ local_vertlevel_list(i) = i
+ end do
+ else
+ allocate(local_vertlevel_list(0))
+ end if
+ allocate(needed_vertlevel_list(nVertLevels+1))
+ do i=1,nVertLevels+1
+ needed_vertlevel_list(i) = i
+ end do
+ call dmpar_get_owner_list(domain % dminfo, &
+ size(local_vertlevel_list), size(needed_vertlevel_list), &
+ local_vertlevel_list, needed_vertlevel_list, &
+ sendVertLevelP1List, recvVertLevelP1List)
+
+ deallocate(local_vertlevel_list)
+ deallocate(needed_vertlevel_list)
+
+
!
! Read and distribute all fields given ownership lists and exchange lists (maybe already in block?)
!
@@ -764,7 +787,7 @@
readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
readVertLevelStart, nReadVertLevels, &
sendCellList, recvCellList, sendEdgeList, recvEdgeList, sendVertexList, recvVertexList, &
- sendVertLevelList, recvVertLevelList)
+ sendVertLevelList, recvVertLevelList, sendVertLevelP1List, recvVertLevelP1List)
call io_input_finalize(input_obj, domain % dminfo)
@@ -804,7 +827,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnCell % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -812,7 +836,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnCell % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -820,7 +845,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnCell % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
end if
end do
@@ -834,7 +860,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnEdge % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -842,7 +869,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnEdge % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
end if
end do
@@ -854,7 +882,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnEdge % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
end if
end do
@@ -868,7 +897,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnVertex % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -876,7 +906,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnVertex % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
end if
end do
@@ -968,7 +999,8 @@
sendCellsList, recvCellsList, &
sendEdgesList, recvEdgesList, &
sendVerticesList, recvVerticesList, &
- sendVertLevelsList, recvVertLevelsList)
+ sendVertLevelsList, recvVertLevelsList, &
+ sendVertLevelsP1List, recvVertLevelsP1List)
implicit none
@@ -981,6 +1013,7 @@
type (exchange_list), pointer :: sendEdgesList, recvEdgesList
type (exchange_list), pointer :: sendVerticesList, recvVerticesList
type (exchange_list), pointer :: sendVertLevelsList, recvVertLevelsList
+ type (exchange_list), pointer :: sendVertLevelsP1List, recvVertLevelsP1List
type (field1dInteger) :: int1d
type (field2dInteger) :: int2d
@@ -1068,6 +1101,33 @@
end subroutine io_input_get_dimension
+ subroutine io_input_field0dReal(input_obj, field)
+
+ implicit none
+
+ type (io_input_object), intent(in) :: input_obj
+ type (field0dReal), intent(inout) :: field
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+ integer :: varID
+ integer, dimension(1) :: start1, count1
+
+ start1(1) = 1
+ count1(1) = 1
+
+#include "input_field0dreal.inc"
+
+#if (RKIND == 8)
+ nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#else
+ nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
+#endif
+
+ end subroutine io_input_field0dReal
+
+
subroutine io_input_field1dReal(input_obj, field)
implicit none
Modified: branches/atmos_physics/src/framework/module_io_output.F
===================================================================
--- branches/atmos_physics/src/framework/module_io_output.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/module_io_output.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -23,10 +23,12 @@
type (exchange_list), pointer :: sendEdgesList, recvEdgesList
type (exchange_list), pointer :: sendVerticesList, recvVerticesList
type (exchange_list), pointer :: sendVertLevelsList, recvVertLevelsList
+ type (exchange_list), pointer :: sendVertLevelsP1List, recvVertLevelsP1List
end type io_output_object
interface io_output_field
+ module procedure io_output_field0dReal
module procedure io_output_field1dReal
module procedure io_output_field2dReal
module procedure io_output_field3dReal
@@ -65,6 +67,8 @@
nullify(output_obj % recvVerticesList)
nullify(output_obj % sendVertLevelsList)
nullify(output_obj % recvVertLevelsList)
+ nullify(output_obj % sendVertLevelsP1List)
+ nullify(output_obj % recvVertLevelsP1List)
output_obj % validExchangeLists = .false.
#include "output_dim_inits.inc"
@@ -109,6 +113,7 @@
integer, dimension(:), pointer :: neededEdgeList
integer, dimension(:), pointer :: neededVertexList
integer, dimension(:), pointer :: neededVertLevelList
+ integer, dimension(:), pointer :: neededVertLevelP1List
integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, verticesOnCell, &
cellsOnEdge, verticesOnEdge, edgesOnEdge, cellsOnVertex, edgesOnVertex
integer, dimension(:,:), pointer :: cellsOnCell_save, edgesOnCell_save, verticesOnCell_save, &
@@ -187,8 +192,12 @@
domain % blocklist % mesh % edgesOnEdge % array(j,i))
end do
do j=domain % blocklist % mesh % nEdgesOnEdge % array(i)+1,2*domain % blocklist % mesh % maxEdges
- edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
+ if(domain % blocklist % mesh % nEdgesOnEdge % array(i) .eq. 0) then
+ edgesOnEdge(j,i) = domain % blocklist % mesh % nEdgesSolve + 1
+ else
+ edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
domain % blocklist % mesh % nEdgesOnEdge % array(i))
+ endif
end do
end do
do i=1,domain % blocklist % mesh % nVerticesSolve
@@ -205,6 +214,7 @@
allocate(neededEdgeList(nEdgesGlobal))
allocate(neededVertexList(nVerticesGlobal))
allocate(neededVertLevelList(nVertLevelsGlobal))
+ allocate(neededVertLevelP1List(nVertLevelsGlobal+1))
do i=1,nCellsGlobal
neededCellList(i) = i
end do
@@ -217,11 +227,15 @@
do i=1,nVertLevelsGlobal
neededVertLevelList(i) = i
end do
+ do i=1,nVertLevelsGlobal+1
+ neededVertLevelP1List(i) = i
+ end do
else
allocate(neededCellList(0))
allocate(neededEdgeList(0))
allocate(neededVertexList(0))
allocate(neededVertLevelList(0))
+ allocate(neededVertLevelP1List(0))
end if
if (.not. output_obj % validExchangeLists) then
@@ -245,6 +259,11 @@
neededVertLevelList, neededVertLevelList, &
output_obj % sendVertLevelsList, output_obj % recvVertLevelsList)
+ call dmpar_get_owner_list(domain % dminfo, &
+ size(neededVertLevelP1List), size(neededVertLevelP1List), &
+ neededVertLevelP1List, neededVertLevelP1List, &
+ output_obj % sendVertLevelsP1List, output_obj % recvVertLevelsP1List)
+
output_obj % validExchangeLists = .true.
end if
@@ -336,6 +355,35 @@
end subroutine io_output_init
+ subroutine io_output_field0dReal(output_obj, field)
+
+ implicit none
+
+ type (io_output_object), intent(in) :: output_obj
+ type (field0dReal), intent(inout) :: field
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+ integer :: varID
+ integer, dimension(1) :: start1, count1
+
+ start1(1) = 1
+ count1(1) = 1
+
+#include "output_field0dreal.inc"
+
+#if (RKIND == 8)
+ nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#else
+ nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
+#endif
+
+ nferr = nf_sync(output_obj % wr_ncid)
+
+ end subroutine io_output_field0dReal
+
+
subroutine io_output_field1dReal(output_obj, field)
implicit none
Modified: branches/atmos_physics/src/framework/module_zoltan_interface.F
===================================================================
--- branches/atmos_physics/src/framework/module_zoltan_interface.F 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/framework/module_zoltan_interface.F 2010-05-20 21:18:40 UTC (rev 295)
@@ -1,9 +1,10 @@
module zoltan_interface
use zoltan
- use mpi
implicit none
+ include 'mpif.h'
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Data for reordering cells
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: branches/atmos_physics/src/operators/Makefile
===================================================================
--- branches/atmos_physics/src/operators/Makefile 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/operators/Makefile 2010-05-20 21:18:40 UTC (rev 295)
@@ -10,10 +10,9 @@
module_vector_reconstruction.o:
clean:
- $(RM) *.o *.mod libops.a
+ $(RM) *.o *.mod *.f90 libops.a
.F.o:
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
- $(RM) $*.f90
Modified: branches/atmos_physics/src/registry/gen_inc.c
===================================================================
--- branches/atmos_physics/src/registry/gen_inc.c 2010-05-20 20:33:43 UTC (rev 294)
+++ branches/atmos_physics/src/registry/gen_inc.c 2010-05-20 21:18:40 UTC (rev 295)
@@ -387,10 +387,20 @@
fortprintf(fd, " allocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(g %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -414,16 +424,28 @@
else {
fortprintf(fd, " allocate(g %% %s)</font>
<font color="black">", var_ptr->name_in_code);
fortprintf(fd, " allocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " allocate(g %% %s %% array(", var_ptr->name_in_code);
- dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
- dimlist_ptr = dimlist_ptr->next;
- while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (var_ptr->ndims > 0) {
+ fortprintf(fd, " allocate(g %% %s %% array(", var_ptr->name_in_code);
+ dimlist_ptr = var_ptr->dimlist;
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
dimlist_ptr = dimlist_ptr->next;
+ while (dimlist_ptr) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
+ dimlist_ptr = dimlist_ptr->next;
+ }
+ fortprintf(fd, "))</font>
<font color="red">");
+
}
- fortprintf(fd, "))</font>
<font color="red">");
-
if (var_ptr->iostreams & INPUT0)
fortprintf(fd, " g %% %s %% ioinfo %% input = .true.</font>
<font color="gray">", var_ptr->name_in_code);
else
@@ -473,9 +495,15 @@
fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="red">", var_ptr2->super_array);
}
else {
- fortprintf(fd, " deallocate(g %% %s %% array)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
- fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="blue">", var_ptr->name_in_code);
+ if (var_ptr->ndims > 0) {
+ fortprintf(fd, " deallocate(g %% %s %% array)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="blue">", var_ptr->name_in_code);
+ }
+ else {
+ fortprintf(fd, " deallocate(g %% %s %% ioinfo)</font>
<font color="blue">", var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(g %% %s)</font>
<font color="black"></font>
<font color="gray">", var_ptr->name_in_code);
+ }
var_ptr = var_ptr->next;
}
}
@@ -508,10 +536,20 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(s %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -538,10 +576,30 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
fortprintf(fd, " allocate(s %% %s %% array(", var_ptr->name_in_code);
dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->constant_value < 0) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
+ }
+ else {
+ fortprintf(fd, "%i", dimlist_ptr->dim->constant_value);
+ }
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->constant_value < 0) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !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);
+ }
+ else {
+ fortprintf(fd, ", %i", dimlist_ptr->dim->constant_value);
+ }
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -976,10 +1034,7 @@
fortprintf(fd, ", block %% mesh %% %s, &</font>
<font color="red">", lastdim->dim->name_in_code);
if (is_derived_dim(lastdim->dim->name_in_code)) {
- split_derived_dim_string(lastdim->dim->name_in_code, &cp1, &cp2);
- fortprintf(fd, " send%sList, recv%sList)</font>
<font color="red">", cp1, cp1);
- free(cp1);
- free(cp2);
+ fortprintf(fd, " send%sList, recv%sList)</font>
<font color="black">", lastdim->dim->name_in_file+1, lastdim->dim->name_in_file+1);
}
else
fortprintf(fd, " send%sList, recv%sList)</font>
<font color="gray">", lastdim->dim->name_in_code+1, lastdim->dim->name_in_code+1);
@@ -1026,7 +1081,10 @@
fortprintf(fd, " deallocate(super_%s%id)</font>
<font color="red">", vtype, var_ptr->ndims);
}
else {
- fortprintf(fd, " block %% time_levs(1) %% state %% %s %% scalar = %s%id %% scalar</font>
<font color="blue">", var_ptr->name_in_code, vtype, var_ptr->ndims);
+ if (var_ptr->timedim)
+ fortprintf(fd, " block %% time_levs(1) %% state %% %s %% scalar = %s%id %% scalar</font>
<font color="blue">", var_ptr->name_in_code, vtype, var_ptr->ndims);
+ else
+ fortprintf(fd, " block %% mesh %% %s %% scalar = %s%id %% scalar</font>
<font color="black">", var_ptr->name_in_code, vtype, var_ptr->ndims);
}
fortprintf(fd, " end if</font>
<font color="black"></font>
<font color="gray">");
@@ -1086,10 +1144,10 @@
/*
- * Generate code to read 1d, 2d, 3d time-invariant fields
+ * Generate code to read 0d, 1d, 2d, 3d time-invariant fields
*/
for(j=0; j<2; j++) {
- for(i=1; i<=3; i++) {
+ for(i=0; i<=3; i++) {
if (j == 0) {
sprintf(fname, "input_field%idinteger.inc", i);
ivtype = INTEGER;
@@ -1486,7 +1544,7 @@
if (is_derived_dim(lastdim->dim->name_in_code)) {
split_derived_dim_string(lastdim->dim->name_in_code, &cp1, &cp2);
fortprintf(fd, ", n%sGlobal%s, &</font>
<font color="red">", cp1, cp2);
- fortprintf(fd, " output_obj %% send%sList, output_obj %% recv%sList)</font>
<font color="blue">", cp1, cp1);
+ fortprintf(fd, " output_obj %% send%sList, output_obj %% recv%sList)</font>
<font color="gray">", lastdim->dim->name_in_file+1, lastdim->dim->name_in_file+1);
free(cp1);
free(cp2);
}
@@ -1497,7 +1555,10 @@
}
else {
fortprintf(fd, " %s%id %% ioinfo %% fieldName = \'%s\'</font>
<font color="red">", vtype, var_ptr->ndims, var_ptr->name_in_file);
- fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% time_levs(1) %% state %% %s %% scalar</font>
<font color="blue">", vtype, var_ptr->ndims, var_ptr->name_in_code);
+ if (var_ptr->timedim)
+ fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% time_levs(1) %% state %% %s %% scalar</font>
<font color="blue">", vtype, var_ptr->ndims, var_ptr->name_in_code);
+ else
+ fortprintf(fd, " %s%id %% scalar = domain %% blocklist %% mesh %% %s %% scalar</font>
<font color="gray">", vtype, var_ptr->ndims, var_ptr->name_in_code);
}
if (var_ptr->timedim)
@@ -1518,10 +1579,10 @@
/*
- * Generate code to write 1d, 2d, 3d time-invariant fields
+ * Generate code to write 0d, 1d, 2d, 3d time-invariant fields
*/
for(j=0; j<2; j++) {
- for(i=1; i<=3; i++) {
+ for(i=0; i<=3; i++) {
if (j == 0) {
sprintf(fname, "output_field%idinteger.inc", i);
ivtype = INTEGER;
</font>
</pre>