<p><b>mpetersen@lanl.gov</b> 2010-03-16 09:10:19 -0600 (Tue, 16 Mar 2010)</p><p>Updates to core_ocean so ocean model runs in new directory structure. Also cleaned up initial conditions in module_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- trunk/mpas/src/core_ocean/module_test_cases.F 2010-03-12 18:22:27 UTC (rev 136)
+++ trunk/mpas/src/core_ocean/module_test_cases.F 2010-03-16 15:10:19 UTC (rev 137)
@@ -21,8 +21,12 @@
type (domain_type), intent(inout) :: domain
- integer :: i
+ integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src
+ real (kind=RKIND), dimension(:), pointer :: rho, alpha
+ 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'
@@ -91,6 +95,65 @@
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
+
+ u_src => block_ptr % mesh % u_src % array
+ rho => block_ptr % mesh % rho % array
+ alpha => block_ptr % mesh % alpha % array
+
+ nCells = block_ptr % mesh % nCells
+ nEdges = block_ptr % mesh % nEdges
+ nVertices = block_ptr % mesh % nVertices
+ nVertLevels = block_ptr % mesh % nVertLevels
+
+
+ ! Momentum forcing u_src
+ ! for shallow water test cases:
+ ! u_src = 0.0
+ ! for rectangular basin:
+ do iEdge=1,nEdges
+ u_src(1:nVertLevels, iEdge) = u_src(1, iEdge) / nVertLevels
+ end do
+
+ ! define the density for multiple layers
+ delta_rho=0.1
+ do iLevel = 1,nVertLevels
+ rho(iLevel) = delta_rho*(iLevel-1)
+ enddo
+ rho = rho - sum(rho)/nVertLevels + 1000.0
+ alpha = 1/rho
+
+ 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), &
+ 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
+ end do
+ ! mrp 100121 end
+
end subroutine setup_sw_test_case
@@ -114,7 +177,7 @@
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
real (kind=RKIND), parameter :: alpha = pii/4.0
- integer :: iCell, iEdge, iVtx, iLevel !mrp on last one
+ integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: r, u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -151,9 +214,6 @@
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
-#ifdef EXPAND_LEVELS
- state % u % array(2:EXPAND_LEVELS, iEdge) = state % u % array(1,iEdge)
-#endif
end do
deallocate(psiVertex)
@@ -167,33 +227,8 @@
else
state % h % array(1,iCell) = 0.0
end if
-#ifdef EXPAND_LEVELS
- ! mrp 100125:
- !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
- ! make the total thickness equal to the single-layer thickness:
- state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
- / grid % nVertLevels
- ! mrp 100125 end
-#else
- ! mrp 100125:
- ! make the total thickness equal to the single-layer thickness:
- !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
- ! mrp 100125 end
-#endif
end do
- ! mrp 100121:
- ! define the density
- do iLevel = 1,grid % nVertLevels
- grid % rho % array(iLevel) = (iLevel-1)*10
- enddo
- grid % rho % array = grid % rho % array &
- - sum(grid % rho % array)/grid % nVertLevels + 1000.0
-print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
- grid % alpha % array = 1/grid % rho % array
- grid % u_src % array = 0.0
- ! mrp 100121 end
-
end subroutine sw_test_case_1
@@ -217,7 +252,6 @@
real (kind=RKIND), parameter :: alpha = 0.0
integer :: iCell, iEdge, iVtx
- integer :: iLevel !mrp
real (kind=RKIND) :: u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -256,9 +290,6 @@
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
-#ifdef EXPAND_LEVELS
- state % u % array(2:EXPAND_LEVELS, iEdge) = state % u % array(1,iEdge)
-#endif
end do
deallocate(psiVertex)
@@ -288,33 +319,8 @@
)**2.0 &
) / &
gravity
-#ifdef EXPAND_LEVELS
- ! mrp 100125:
- !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
- ! make the total thickness equal to the single-layer thickness:
- state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
- / grid % nVertLevels
- ! mrp 100125 end
-#else
- ! mrp 100125:
- ! make the total thickness equal to the single-layer thickness:
- !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
- ! mrp 100125 end
-#endif
end do
- ! mrp 100121:
- ! define the density
- do iLevel = 1,grid % nVertLevels
- grid % rho % array(iLevel) = (iLevel-1)*10
- enddo
- grid % rho % array = grid % rho % array &
- - sum(grid % rho % array)/grid % nVertLevels + 1000.0
-print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
- grid % alpha % array = 1/grid % rho % array
- grid % u_src % array = 0.0
- ! mrp 100121 end
-
end subroutine sw_test_case_2
@@ -341,7 +347,7 @@
real (kind=RKIND), parameter :: rr = pii/9.0
real (kind=RKIND), parameter :: alpha = 0.0
- integer :: iCell, iEdge, iVtx, iLevel
+ integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: r, u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -379,9 +385,6 @@
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
-#ifdef EXPAND_LEVELS
- state % u % array(2:EXPAND_LEVELS, iEdge) = state % u % array(1,iEdge)
-#endif
end do
deallocate(psiVertex)
@@ -438,33 +441,8 @@
) / &
gravity
state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
-#ifdef EXPAND_LEVELS
- ! mrp 100125:
- !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
- ! make the total thickness equal to the single-layer thickness:
- state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
- / grid % nVertLevels
- ! mrp 100125 end
-#else
- ! mrp 100125:
- ! make the total thickness equal to the single-layer thickness:
- !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
- ! mrp 100125 end
-#endif
end do
- ! mrp 100121:
- ! define the density
- do iLevel = 1,grid % nVertLevels
- grid % rho % array(iLevel) = 0.0 + (iLevel-1)*10
- enddo
- grid % rho % array = grid % rho % array &
- - sum(grid % rho % array)/grid % nVertLevels + 1000.0
-print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
- grid % alpha % array = 1/grid % rho % array
- grid % u_src % array = 0.0
- ! mrp 100121 end
-
end subroutine sw_test_case_5
@@ -487,7 +465,7 @@
real (kind=RKIND), parameter :: K = 7.848e-6
real (kind=RKIND), parameter :: R = 4.0
- integer :: iCell, iEdge, iVtx, iLevel
+ integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -524,9 +502,6 @@
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
-#ifdef EXPAND_LEVELS
- state % u % array(2:EXPAND_LEVELS, iEdge) = state % u % array(1,iEdge)
-#endif
end do
deallocate(psiVertex)
@@ -538,33 +513,8 @@
a*a*BB(grid%latCell%array(iCell)) * cos(R*grid%lonCell%array(iCell)) + &
a*a*CC(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
) / gravity
-#ifdef EXPAND_LEVELS
- ! mrp 100125:
- !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
- ! make the total thickness equal to the single-layer thickness:
- state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
- / grid % nVertLevels
- ! mrp 100125 end
-#else
- ! mrp 100125:
- ! make the total thickness equal to the single-layer thickness:
- !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
- ! mrp 100125 end
-#endif
end do
- ! mrp 100121:
- ! define the density
- do iLevel = 1,grid % nVertLevels
- grid % rho % array(iLevel) = (iLevel-1)*10
- enddo
- grid % rho % array = grid % rho % array &
- - sum(grid % rho % array)/grid % nVertLevels + 1000.0
-print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
- grid % alpha % array = 1/grid % rho % array
- grid % u_src % array = 0.0
- ! mrp 100121 end
-
end subroutine sw_test_case_6
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2010-03-12 18:22:27 UTC (rev 136)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2010-03-16 15:10:19 UTC (rev 137)
@@ -26,6 +26,7 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ integer errorcode,ierr
type (block_type), pointer :: block
@@ -40,6 +41,13 @@
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
@@ -125,6 +133,7 @@
block => domain % blocklist
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)
Modified: trunk/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_interface.F 2010-03-12 18:22:27 UTC (rev 136)
+++ trunk/mpas/src/core_ocean/mpas_interface.F 2010-03-16 15:10:19 UTC (rev 137)
@@ -30,8 +30,12 @@
! call init_reconstruct(block_ptr % mesh)
! xsad 10-02-09 end
+! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
+! input arguement into mpas_init. Ask about that later. For now, there will be
+! no initial statistics write.
+
! call timer_start("global diagnostics")
-! call computeGlobalDiagnostics(domain % dminfo, block_ptr % time_levs(1) % state, block_ptr % mesh, 0, dt)
+! call computeGlobalDiagnostics(domain % dminfo, block % time_levs(1) % state, mesh, 0, dt)
! call timer_stop("global diagnostics")
! ! xsad 10-02-08 end
! call output_state_init(output_obj, domain, "OUTPUT")
@@ -48,7 +52,7 @@
integer, intent(out) :: ivalue
if (index(key,'STORAGE_FACTOR') /= 0) then
- ivalue = 1
+ ivalue = 2
end if
end subroutine mpas_query
@@ -58,33 +62,37 @@
use grid_types
use time_integration
+ use timer
+ use global_diagnostics
implicit none
type (domain_type), intent(inout) :: domain
integer, intent(in) :: itimestep
real (kind=RKIND), intent(in) :: dt
+ type (block_type), pointer :: block_ptr
call timestep(domain, dt)
-! ! mrp 100120:
-! if (mod(itimestep, config_stats_interval) == 0) then
-! ! xsad 10-02-08:
-! !call write_stats(domain, itimestep, dt)
-! block_ptr => domain % blocklist
-! if(associated(block_ptr % next)) then
-! write(0,*) 'Error: computeGlobalDiagnostics assumes that there is only one block per processor.'
-! call dmpar_abort(domain % dminfo)
-! end if
-!
-! call timer_start("global diagnostics")
-! call computeGlobalDiagnostics(domain % dminfo, block_ptr % time_levs(2) % state, block_ptr % mesh, &
-! itimestep, dt)
-! call timer_stop("global diagnostics")
-! ! xsad 10-02-08 end
-! end if
-! ! mrp 100120 end
+ ! mrp 100120:
+ if (mod(itimestep, config_stats_interval) == 0) then
+ ! xsad 10-02-08:
+ !call write_stats(domain, itimestep, dt)
+ block_ptr => domain % blocklist
+ if(associated(block_ptr % next)) then
+ write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+ 'that there is only one block per processor.'
+ end if
+ call timer_start("global diagnostics")
+ call computeGlobalDiagnostics(domain % dminfo, &
+ block_ptr % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global diagnostics")
+ ! xsad 10-02-08 end
+ end if
+ ! mrp 100120 end
+
end subroutine mpas_timestep
</font>
</pre>