<p><b>mpetersen@lanl.gov</b> 2010-11-15 16:08:12 -0700 (Mon, 15 Nov 2010)</p><p>BRANCH COMMIT Merge trunk to topography branch, notably the new driver structure.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -105,6 +105,18 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+gfortran-serial:
+ ( make all \
+ "FC = gfortran" \
+ "CC = gcc" \
+ "SFC = gfortran" \
+ "SCC = gcc" \
+ "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
+ "CFLAGS = -O3 -m64" \
+ "LDFLAGS = -O3 -m64" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
g95:
( make all \
"FC = mpif90" \
Modified: branches/ocean_projects/topography2_mrp/namelist.input.sw
===================================================================
--- branches/ocean_projects/topography2_mrp/namelist.input.sw 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/namelist.input.sw 2010-11-15 23:08:12 UTC (rev 617)
@@ -5,7 +5,10 @@
config_ntimesteps = 7500
config_output_interval = 500
config_stats_interval = 0
- config_visc = 0.0
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
config_tracer_adv_order = 2
config_positive_definite = .false.
Modified: branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,9 +1,9 @@
.SUFFIXES: .F .o
-OBJS = module_test_cases.o \
+OBJS = module_mpas_core.o \
+ module_test_cases.o \
module_time_integration.o \
- module_advection.o \
- mpas_interface.o
+ module_advection.o
all: core_hyd
@@ -16,7 +16,7 @@
module_advection.o:
-mpas_interface.o: module_advection.o module_test_cases.o module_time_integration.o
+module_mpas_core.o: module_advection.o module_test_cases.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Registry 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/Registry 2010-11-15 23:08:12 UTC (rev 617)
@@ -18,6 +18,7 @@
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic true
namelist integer sw_model config_mp_physics 0
+namelist real sw_model config_apvm_upwinding 0.5
namelist integer dimensions config_nvertlevels 26
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
Copied: branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_mpas_core.F (from rev 616, trunk/mpas/src/core_hyd_atmos/module_mpas_core.F)
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_mpas_core.F (rev 0)
+++ branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_mpas_core.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -0,0 +1,201 @@
+module mpas_core
+
+ use mpas_framework
+
+ type (io_output_object) :: restart_obj
+ integer :: restart_frame
+
+
+ contains
+
+
+ subroutine mpas_core_init(domain)
+
+ use configure
+ use grid_types
+ use test_cases
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block
+
+
+ if (.not. config_do_restart) call setup_hyd_test_case(domain)
+
+ !
+ ! Initialize core
+ !
+ dt = config_dt
+ block => domain % blocklist
+ do while (associated(block))
+ call mpas_init_block(block, block % mesh, dt)
+ block => block % next
+ end do
+
+ restart_frame = 1
+
+ end subroutine mpas_core_init
+
+
+ subroutine mpas_init_block(block, mesh, dt)
+
+ use grid_types
+ use advection
+ use time_integration
+ use RBF_interpolation
+ use vector_reconstruction
+
+ implicit none
+
+ type (block_type), intent(inout) :: block
+ type (mesh_type), intent(inout) :: mesh
+ real (kind=RKIND), intent(in) :: dt
+
+
+ call compute_solver_constants(block % state % time_levs(1) % state, mesh)
+ call compute_state_diagnostics(block % state % time_levs(1) % state, mesh)
+ call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+ call initialize_advection_rk(mesh)
+ call rbfInterp_initialize(mesh)
+ call init_reconstruct(mesh)
+ call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
+
+ end subroutine mpas_init_block
+
+
+ subroutine mpas_core_run(domain, output_obj, output_frame)
+
+ use grid_types
+ use io_output
+ use timer
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+
+ integer :: ntimesteps, itimestep
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block_ptr
+
+ ! Eventually, dt should be domain specific
+ dt = config_dt
+ ntimesteps = config_ntimesteps
+
+ call write_output_frame(output_obj, output_frame, domain)
+
+ ! During integration, time level 1 stores the model state at the beginning of the
+ ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
+ do itimestep = 1,ntimesteps
+ write(0,*) 'Doing timestep ', itimestep
+ call timer_start("time integration")
+ call mpas_timestep(domain, itimestep, dt)
+ call timer_stop("time integration")
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels_state(domain % blocklist % state)
+
+ if (mod(itimestep, config_output_interval) == 0) then
+ call write_output_frame(output_obj, output_frame, domain)
+ end if
+ if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
+ call output_state_for_domain(restart_obj, domain, restart_frame)
+ restart_frame = restart_frame + 1
+ end if
+ end do
+
+ end subroutine mpas_core_run
+
+
+ subroutine write_output_frame(output_obj, output_frame, domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain and write model state to output file
+ !
+ ! Input/Output: domain - contains model state; diagnostic field are computed
+ ! before returning
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+ use io_output
+
+ implicit none
+
+ integer, intent(inout) :: output_frame
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+
+ integer :: i, j, k
+ integer :: eoe
+ type (block_type), pointer :: block_ptr
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
+ block_ptr => block_ptr % next
+ end do
+
+ call output_state_for_domain(output_obj, domain, output_frame)
+ output_frame = output_frame + 1
+
+ end subroutine write_output_frame
+
+
+ subroutine compute_output_diagnostics(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain
+ !
+ ! Input: state - contains model prognostic fields
+ ! grid - contains grid metadata
+ !
+ ! Output: state - upon returning, diagnostic fields will have be computed
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, eoe
+ integer :: iEdge, k
+
+ end subroutine compute_output_diagnostics
+
+
+ subroutine mpas_timestep(domain, itimestep, dt)
+
+ use grid_types
+ use time_integration
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: dt
+
+ call timestep(domain, dt)
+
+ end subroutine mpas_timestep
+
+
+ subroutine mpas_core_finalize(domain)
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
+
+ end subroutine mpas_core_finalize
+
+end module mpas_core
Modified: branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_time_integration.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/module_time_integration.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1988,7 +1988,7 @@
!
do iEdge = 1,nEdges
do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
end do
end do
@@ -2027,7 +2027,7 @@
!
do iEdge = 1,nEdges
do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
end do
end do
Deleted: branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/mpas_interface.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_hyd_atmos/mpas_interface.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,74 +0,0 @@
-subroutine mpas_setup_test_case(domain)
-
- use grid_types
- use test_cases
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- call setup_hyd_test_case(domain)
-
-end subroutine mpas_setup_test_case
-
-
-subroutine mpas_init(block, mesh, dt)
-
- use grid_types
- use advection
- use time_integration
- use RBF_interpolation
- use vector_reconstruction
-
- implicit none
-
- type (block_type), intent(inout) :: block
- type (mesh_type), intent(inout) :: mesh
- real (kind=RKIND), intent(in) :: dt
-
- call compute_solver_constants(block % state % time_levs(1) % state, mesh)
- call compute_state_diagnostics(block % state % time_levs(1) % state, mesh)
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
- call initialize_advection_rk(mesh)
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(block % state % time_levs(1) % state, mesh)
-
-end subroutine mpas_init
-
-
-subroutine mpas_query(key, ivalue)
-
- implicit none
-
- character (len=256), intent(in) :: key
- integer, intent(out) :: ivalue
-
- if (index(key,'STORAGE_FACTOR') /= 0) then
- ivalue = 1
- end if
-
-end subroutine mpas_query
-
-
-subroutine mpas_timestep(domain, itimestep, dt)
-
- use grid_types
- use time_integration
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- integer, intent(in) :: itimestep
- real (kind=RKIND), intent(in) :: dt
-
- call timestep(domain, dt)
-
-end subroutine mpas_timestep
-
-
-subroutine mpas_finalize()
-
- implicit none
-
-end subroutine mpas_finalize
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,10 +1,10 @@
.SUFFIXES: .F .o
-OBJS = module_test_cases.o \
+OBJS = module_mpas_core.o \
+ module_test_cases.o \
module_advection.o \
module_time_integration.o \
- module_global_diagnostics.o \
- mpas_interface.o
+ module_global_diagnostics.o
all: core_hyd
@@ -19,7 +19,7 @@
module_global_diagnostics.o:
-mpas_interface.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_time_integration.o
+module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-15 23:08:12 UTC (rev 617)
@@ -19,6 +19,7 @@
namelist real hmix config_h_mom_eddy_visc4 0.0
namelist real hmix config_h_tracer_eddy_diff2 0.0
namelist real hmix config_h_tracer_eddy_diff4 0.0
+namelist real hmix config_apvm_upwinding 0.5
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
namelist real vmix config_vert_viscosity 2.5e-4
Copied: branches/ocean_projects/topography2_mrp/src/core_ocean/module_mpas_core.F (from rev 616, trunk/mpas/src/core_ocean/module_mpas_core.F)
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_mpas_core.F (rev 0)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_mpas_core.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -0,0 +1,414 @@
+module mpas_core
+
+ use mpas_framework
+ use dmpar
+ use test_cases
+
+ type (io_output_object) :: restart_obj
+ integer :: restart_frame
+
+
+ contains
+
+
+ subroutine mpas_core_init(domain)
+
+ use configure
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block
+ type (dm_info) :: dminfo
+
+
+ if (.not. config_do_restart) call setup_sw_test_case(domain)
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+ print *, ' Using isopycnal coordinates'
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ print *, ' Using z-level coordinates'
+ call init_ZLevel(domain)
+ else
+ print *, ' Incorrect choice of config_vert_grid_type:',&
+ config_vert_grid_type
+ call dmpar_abort(dminfo)
+ endif
+
+ call compute_maxLevel(domain)
+
+ !
+ ! Initialize core
+ !
+ dt = config_dt
+ block => domain % blocklist
+ do while (associated(block))
+ call mpas_init_block(block, block % mesh, dt)
+ block => block % next
+ end do
+
+ ! 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 % state % time_levs(1) % state, mesh, 0, dt)
+ ! call timer_stop("global diagnostics")
+ ! call output_state_init(output_obj, domain, "OUTPUT")
+ ! call write_output_frame(output_obj, domain)
+
+ restart_frame = 1
+
+ end subroutine mpas_core_init
+
+
+ subroutine mpas_init_block(block, mesh, dt)
+
+ use grid_types
+ use time_integration
+ use RBF_interpolation
+ use vector_reconstruction
+
+ implicit none
+
+ type (block_type), intent(inout) :: block
+ type (mesh_type), intent(inout) :: mesh
+ real (kind=RKIND), intent(in) :: dt
+
+
+ call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+
+ call rbfInterp_initialize(mesh)
+ call init_reconstruct(mesh)
+ call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
+
+ end subroutine mpas_init_block
+
+
+ subroutine mpas_core_run(domain, output_obj, output_frame)
+
+ use grid_types
+ use io_output
+ use timer
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+
+ integer :: ntimesteps, itimestep
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block_ptr
+
+ ! Eventually, dt should be domain specific
+ dt = config_dt
+ ntimesteps = config_ntimesteps
+
+ call write_output_frame(output_obj, output_frame, domain)
+
+ ! During integration, time level 1 stores the model state at the beginning of the
+ ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
+ do itimestep = 1,ntimesteps
+ write(0,*) 'Doing timestep ', itimestep
+ call timer_start("time integration")
+ call mpas_timestep(domain, itimestep, dt)
+ call timer_stop("time integration")
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels_state(domain % blocklist % state)
+
+ if (mod(itimestep, config_output_interval) == 0) then
+ call write_output_frame(output_obj, output_frame, domain)
+ end if
+ if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
+ call output_state_for_domain(restart_obj, domain, restart_frame)
+ restart_frame = restart_frame + 1
+ end if
+ end do
+
+ end subroutine mpas_core_run
+
+
+ subroutine write_output_frame(output_obj, output_frame, domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain and write model state to output file
+ !
+ ! Input/Output: domain - contains model state; diagnostic field are computed
+ ! before returning
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+ use io_output
+
+ implicit none
+
+ integer, intent(inout) :: output_frame
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+
+ integer :: i, j, k
+ integer :: eoe
+ type (block_type), pointer :: block_ptr
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
+ block_ptr => block_ptr % next
+ end do
+
+ call output_state_for_domain(output_obj, domain, output_frame)
+ output_frame = output_frame + 1
+
+ end subroutine write_output_frame
+
+
+ subroutine compute_output_diagnostics(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain
+ !
+ ! Input: state - contains model prognostic fields
+ ! grid - contains grid metadata
+ !
+ ! Output: state - upon returning, diagnostic fields will have be computed
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, eoe
+ integer :: iEdge, k
+
+ end subroutine compute_output_diagnostics
+
+
+ subroutine mpas_timestep(domain, itimestep, dt)
+
+ 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)
+
+ if (mod(itimestep, config_stats_interval) == 0) then
+ 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 % state % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global diagnostics")
+ end if
+
+ end subroutine mpas_timestep
+
+
+subroutine init_ZLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell
+ real (kind=RKIND), dimension(:), pointer :: &
+ hZLevel, zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ hZLevel => block % mesh % hZLevel % array
+ zMidZLevel => block % mesh % zMidZLevel % array
+ zTopZLevel => block % mesh % zTopZLevel % array
+ nVertLevels = block % mesh % nVertLevels
+
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ ! Upon restart, the correct hZLevel should be in restart.nc
+ if (.not. config_do_restart) hZLevel = h(:,1)
+
+ ! hZLevel should be in the grid.nc and restart.nc file,
+ ! and h for k=1 must be specified there as well.
+
+ zTopZLevel(1) = 0.0
+ do k = 1,nVertLevels
+ zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
+ zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
+ end do
+
+ block => block % next
+
+ end do
+
+end subroutine init_ZLevel
+
+
+subroutine compute_maxLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
+ real (kind=RKIND) :: centerx, centery
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+ integer, dimension(:), pointer :: &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell, &
+ boundaryVertex, verticesOnEdge
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+ maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+ maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ verticesOnEdge => block % mesh % verticesOnEdge % array
+ boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryCell => block % mesh % boundaryCell % array
+ boundaryVertex => block % mesh % boundaryVertex % array
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+ vertexDegree = block % mesh % vertexDegree
+
+ ! for z-grids, maxLevelCell should be in input state
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
+ maxLevelCell(1:nCells) = nVertLevels
+ endif
+ maxLevelCell(nCells+1) = 0
+
+ ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeTop(nEdges+1) = 0
+
+ ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeBot(iEdge) = &
+ max( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeBot(nEdges+1) = 0
+
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexBot(nVertices+1) = 0
+
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
+
+ ! set boundary edge
+ boundaryEdge=1
+ do iEdge=1,nEdges
+ boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
+ end do
+
+ !
+ ! Find cells and vertices that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
+ boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
+ endif
+ end do
+ end do
+
+ block => block % next
+ end do
+
+ ! Note: We do not update halos on maxLevel* variables. I want the
+ ! outside edge of a halo to be zero on each processor.
+
+end subroutine compute_maxLevel
+
+
+ subroutine mpas_core_finalize(domain)
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
+
+ end subroutine mpas_core_finalize
+
+end module mpas_core
Deleted: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,301 +0,0 @@
-subroutine mpas_setup_test_case(domain)
-
- use grid_types
- use test_cases
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- call setup_sw_test_case(domain)
-
-end subroutine mpas_setup_test_case
-
-
-subroutine mpas_init_domain(domain)
-! Initialize grid variables that are computed from the netcdf input file.
-
- use grid_types
- use dmpar
- use configure
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- type (dm_info) :: dminfo
-
- if (config_vert_grid_type.eq.'isopycnal') then
- print *, ' Using isopycnal coordinates'
- elseif (config_vert_grid_type.eq.'zlevel') then
- print *, ' Using z-level coordinates'
- call init_ZLevel(domain)
- else
- print *, ' Incorrect choice of config_vert_grid_type:',&
- config_vert_grid_type
- call dmpar_abort(dminfo)
- endif
-
-
- call compute_maxLevel(domain)
-
-end subroutine mpas_init_domain
-
-
-subroutine init_ZLevel(domain)
-! Initialize maxLevel and bouncary grid variables.
-
- use grid_types
- use configure
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i, iCell, iEdge, iVertex, k
- type (block_type), pointer :: block
-
- integer :: iTracer, cell
- real (kind=RKIND), dimension(:), pointer :: &
- hZLevel, zMidZLevel, zTopZLevel
- real (kind=RKIND), dimension(:,:), pointer :: h
- integer :: nVertLevels
-
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
-
- h => block % state % time_levs(1) % state % h % array
- hZLevel => block % mesh % hZLevel % array
- zMidZLevel => block % mesh % zMidZLevel % array
- zTopZLevel => block % mesh % zTopZLevel % array
- nVertLevels = block % mesh % nVertLevels
-
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- ! Upon restart, the correct hZLevel should be in restart.nc
- if (.not. config_do_restart) hZLevel = h(:,1)
-
- ! hZLevel should be in the grid.nc and restart.nc file,
- ! and h for k=1 must be specified there as well.
-
- zTopZLevel(1) = 0.0
- do k = 1,nVertLevels
- zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
- zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
- end do
-
- block => block % next
-
- end do
-
-end subroutine init_ZLevel
-
-
-subroutine compute_maxLevel(domain)
-! Initialize maxLevel and bouncary grid variables.
-
- use grid_types
- use configure
- use constants
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i, iCell, iEdge, iVertex, k
- type (block_type), pointer :: block
-
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
-
- integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell, &
- boundaryVertex, verticesOnEdge
-
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
-
- maxLevelCell => block % mesh % maxLevelCell % array
- maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
- maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
- maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
- maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
- cellsOnEdge => block % mesh % cellsOnEdge % array
- cellsOnVertex => block % mesh % cellsOnVertex % array
- verticesOnEdge => block % mesh % verticesOnEdge % array
- boundaryEdge => block % mesh % boundaryEdge % array
- boundaryCell => block % mesh % boundaryCell % array
- boundaryVertex => block % mesh % boundaryVertex % array
-
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertices = block % mesh % nVertices
- nVertLevels = block % mesh % nVertLevels
- vertexDegree = block % mesh % vertexDegree
-
- ! for z-grids, maxLevelCell should be in input state
- ! Isopycnal grid uses all vertical cells
- if (config_vert_grid_type.eq.'isopycnal') then
- maxLevelCell(1:nCells) = nVertLevels
- endif
- maxLevelCell(nCells+1) = 0
-
- ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
- do iEdge=1,nEdges
- maxLevelEdgeTop(iEdge) = &
- min( maxLevelCell(cellsOnEdge(1,iEdge)), &
- maxLevelCell(cellsOnEdge(2,iEdge)) )
- end do
- maxLevelEdgeTop(nEdges+1) = 0
-
- ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
- do iEdge=1,nEdges
- maxLevelEdgeBot(iEdge) = &
- max( maxLevelCell(cellsOnEdge(1,iEdge)), &
- maxLevelCell(cellsOnEdge(2,iEdge)) )
- end do
- maxLevelEdgeBot(nEdges+1) = 0
-
- ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
- do iVertex = 1,nVertices
- maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
- do i=2,vertexDegree
- maxLevelVertexBot(iVertex) = &
- max( maxLevelVertexBot(iVertex), &
- maxLevelCell(cellsOnVertex(i,iVertex)))
- end do
- end do
- maxLevelVertexBot(nVertices+1) = 0
-
- ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
- do iVertex = 1,nVertices
- maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
- do i=2,vertexDegree
- maxLevelVertexTop(iVertex) = &
- min( maxLevelVertexTop(iVertex), &
- maxLevelCell(cellsOnVertex(i,iVertex)))
- end do
- end do
- maxLevelVertexTop(nVertices+1) = 0
-
- ! set boundary edge
- boundaryEdge=1
- do iEdge=1,nEdges
- boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
- end do
-
- !
- ! Find cells and vertices that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if (boundaryEdge(k,iEdge).eq.1) then
- boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
- boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
- boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
- boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
- endif
- end do
- end do
-
- block => block % next
- end do
-
- ! Note: We do not update halos on maxLevel* variables. I want the
- ! outside edge of a halo to be zero on each processor.
-
-end subroutine compute_maxLevel
-
-
-subroutine mpas_init(block, mesh, dt)
-
- use grid_types
- use time_integration
- use RBF_interpolation
- use vector_reconstruction
-
- implicit none
-
- type (block_type), intent(inout) :: block
- type (mesh_type), intent(inout) :: mesh
- real (kind=RKIND), intent(in) :: dt
-
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(block % state % time_levs(1) % state, mesh)
-
-! 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 % state % time_levs(1) % state, mesh, 0, dt)
-! call timer_stop("global diagnostics")
-! call output_state_init(output_obj, domain, "OUTPUT")
-! call write_output_frame(output_obj, domain)
-
-end subroutine mpas_init
-
-
-subroutine mpas_query(key, ivalue)
-
- implicit none
-
- character (len=256), intent(in) :: key
- integer, intent(out) :: ivalue
-
- if (index(key,'STORAGE_FACTOR') /= 0) then
- ivalue = 2
- end if
-
-end subroutine mpas_query
-
-
-subroutine mpas_timestep(domain, itimestep, dt)
-
- 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
-
- call timestep(domain, dt)
-
- if (mod(itimestep, config_stats_interval) == 0) then
- block => domain % blocklist
- if(associated(block % 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 % state % time_levs(2) % state, block % mesh, &
- itimestep, dt)
- call timer_stop("global diagnostics")
- end if
-
-end subroutine mpas_timestep
-
-
-subroutine mpas_finalize()
-
- implicit none
-
-end subroutine mpas_finalize
Modified: branches/ocean_projects/topography2_mrp/src/core_sw/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,10 +1,10 @@
.SUFFIXES: .F .o
-OBJS = module_test_cases.o \
+OBJS = module_mpas_core.o \
+ module_test_cases.o \
module_advection.o \
module_time_integration.o \
- module_global_diagnostics.o \
- mpas_interface.o
+ module_global_diagnostics.o
all: core_sw
@@ -19,7 +19,7 @@
module_global_diagnostics.o:
-mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o
+module_mpas_core.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/ocean_projects/topography2_mrp/src/core_sw/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/Registry 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/Registry 2010-11-15 23:08:12 UTC (rev 617)
@@ -7,11 +7,15 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist integer sw_model config_stats_interval 100
-namelist real sw_model config_visc 0.0
+namelist real sw_model config_h_mom_eddy_visc2 0.0
+namelist real sw_model config_h_mom_eddy_visc4 0.0
+namelist real sw_model config_h_tracer_eddy_diff2 0.0
+namelist real sw_model config_h_tracer_eddy_diff4 0.0
namelist integer sw_model config_thickness_adv_order 2
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
+namelist real sw_model config_apvm_upwinding 0.5
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -140,3 +144,4 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
+
Modified: branches/ocean_projects/topography2_mrp/src/core_sw/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/module_global_diagnostics.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/module_global_diagnostics.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -274,7 +274,7 @@
else
open(fileID, file='GlobalIntegrals.txt',POSITION='append')
endif
- write(fileID,'(1i, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
+ write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
globalKineticEnergy, globalPotentialEnergy
close(fileID)
Copied: branches/ocean_projects/topography2_mrp/src/core_sw/module_mpas_core.F (from rev 616, trunk/mpas/src/core_sw/module_mpas_core.F)
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/module_mpas_core.F (rev 0)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/module_mpas_core.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -0,0 +1,217 @@
+module mpas_core
+
+ use mpas_framework
+
+ type (io_output_object) :: restart_obj
+ integer :: restart_frame
+
+
+ contains
+
+
+ subroutine mpas_core_init(domain)
+
+ use configure
+ use grid_types
+ use test_cases
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block
+
+
+ if (.not. config_do_restart) call setup_sw_test_case(domain)
+
+ !
+ ! Initialize core
+ !
+ dt = config_dt
+ block => domain % blocklist
+ do while (associated(block))
+ call mpas_init_block(block, block % mesh, dt)
+ block => block % next
+ end do
+
+ restart_frame = 1
+
+ end subroutine mpas_core_init
+
+
+ subroutine mpas_init_block(block, mesh, dt)
+
+ use grid_types
+ use time_integration
+ use RBF_interpolation
+ use vector_reconstruction
+
+ implicit none
+
+ type (block_type), intent(inout) :: block
+ type (mesh_type), intent(inout) :: mesh
+ real (kind=RKIND), intent(in) :: dt
+
+
+ call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+
+ call rbfInterp_initialize(mesh)
+ call init_reconstruct(mesh)
+ call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
+
+ end subroutine mpas_init_block
+
+
+ subroutine mpas_core_run(domain, output_obj, output_frame)
+
+ use grid_types
+ use io_output
+ use timer
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+
+ integer :: ntimesteps, itimestep
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block_ptr
+
+ ! Eventually, dt should be domain specific
+ dt = config_dt
+ ntimesteps = config_ntimesteps
+
+ call write_output_frame(output_obj, output_frame, domain)
+
+ ! During integration, time level 1 stores the model state at the beginning of the
+ ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
+ do itimestep = 1,ntimesteps
+ write(0,*) 'Doing timestep ', itimestep
+ call timer_start("time integration")
+ call mpas_timestep(domain, itimestep, dt)
+ call timer_stop("time integration")
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels_state(domain % blocklist % state)
+
+ if (mod(itimestep, config_output_interval) == 0) then
+ call write_output_frame(output_obj, output_frame, domain)
+ end if
+ if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
+ call output_state_for_domain(restart_obj, domain, restart_frame)
+ restart_frame = restart_frame + 1
+ end if
+ end do
+
+ end subroutine mpas_core_run
+
+
+ subroutine write_output_frame(output_obj, output_frame, domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain and write model state to output file
+ !
+ ! Input/Output: domain - contains model state; diagnostic field are computed
+ ! before returning
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+ use io_output
+
+ implicit none
+
+ integer, intent(inout) :: output_frame
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+
+ integer :: i, j, k
+ integer :: eoe
+ type (block_type), pointer :: block_ptr
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
+ block_ptr => block_ptr % next
+ end do
+
+ call output_state_for_domain(output_obj, domain, output_frame)
+ output_frame = output_frame + 1
+
+ end subroutine write_output_frame
+
+
+ subroutine compute_output_diagnostics(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain
+ !
+ ! Input: state - contains model prognostic fields
+ ! grid - contains grid metadata
+ !
+ ! Output: state - upon returning, diagnostic fields will have be computed
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, eoe
+ integer :: iEdge, k
+
+ end subroutine compute_output_diagnostics
+
+
+ subroutine mpas_timestep(domain, itimestep, dt)
+
+ 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)
+
+ if(config_stats_interval .gt. 0) then
+ if(mod(itimestep, config_stats_interval) == 0) then
+ 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 % state % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global_diagnostics")
+ end if
+ end if
+
+ end subroutine mpas_timestep
+
+
+ subroutine mpas_core_finalize(domain)
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
+
+ end subroutine mpas_core_finalize
+
+end module mpas_core
Modified: branches/ocean_projects/topography2_mrp/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/module_time_integration.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/module_time_integration.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -120,6 +120,16 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ end if
+
block => block % next
end do
@@ -253,11 +263,14 @@
circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: u_diffusion, visc
+ real (kind=RKIND) :: r, u_diffusion
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
- visc = config_visc
+
h => s % h % array
u => s % u % array
v => s % v % array
@@ -334,15 +347,6 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
-
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
- ! only valid for visc == constant
- !
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = visc * u_diffusion
-
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
@@ -352,12 +356,13 @@
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)
end do
end do
+
+
#endif
#ifdef NCAR_FORMULATION
@@ -386,6 +391,119 @@
end do
#endif
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for visc == constant
+ if (config_h_mom_eddy_visc2 > 0.0) then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+ end if
+
+ !
+ ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="blue">abla^4 u
+ ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ ! applied recursively.
+ ! strictly only valid for h_mom_eddy_visc4 == constant
+ !
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
+
+ delsq_u(:,:) = 0.0
+
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ end do
+ end do
+
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end do
+ do iVertex=1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k=1,nVertLevels
+ delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ end do
+ end do
+
+ ! Divergence using </font>
<font color="blue">abla^2 u
+ delsq_divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ 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)
+ do k = 1,nVertLevels
+ delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ end do
+ end do
+
+ ! Compute - \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ u_diffusion = ( delsq_divergence(k,cell2) &
+ - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( delsq_vorticity(k,vertex2) &
+ - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = config_h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
+ end if
+
end subroutine compute_tend
@@ -405,7 +523,12 @@
type (mesh_type), intent(in) :: grid
integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
+
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
@@ -422,6 +545,7 @@
tracers => s % tracers % array
cellsOnEdge => grid % cellsOnEdge % array
boundaryCell=> grid % boundaryCell % array
+ boundaryEdge=> grid % boundaryEdge % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
@@ -429,6 +553,7 @@
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
tracer_tend(:,:,:) = 0.0
if (config_tracer_adv_order == 2) then
@@ -556,6 +681,108 @@
endif ! if (config_tracer_adv_order == 2 )
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="blue">abla \phi)
+ !
+ if ( config_h_tracer_eddy_diff2 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1.0/areaCell(cell1)
+ invAreaCell2 = 1.0/areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+ tracer_turb_flux = config_h_tracer_eddy_diff2 &
+ *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
+
+ ! div(h \kappa_2 </font>
<font color="blue">abla \phi) at cell center
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
+ end do
+ end do
+
+ end do
+
+ deallocate(boundaryMask)
+
+ end if
+
+ !
+ ! tracer tendency: del4 horizontal tracer diffusion, &
+ ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="blue">abla \phi)])
+ !
+ if ( config_h_tracer_eddy_diff4 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ ! first del2: div(h </font>
<font color="blue">abla \phi) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
+ - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ end do
+ end do
+
+ end do
+
+ do iCell = 1, grid % nCells
+ r = 1.0 / grid % areaCell % array(iCell)
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
+ end do
+ end do
+
+ ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
+ invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
+ flux = dvEdge(iEdge) * tracer_turb_flux
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ end do
+ enddo
+
+ end do
+
+ deallocate(delsq_tracer)
+ deallocate(boundaryMask)
+
+ end if
+
end subroutine compute_scalar_tend
@@ -779,16 +1006,10 @@
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- 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) <= nVertices) 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
@@ -845,11 +1066,9 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe <= nEdges) 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
@@ -902,26 +1121,25 @@
!
! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
+ ! ( this computes pv_edge at all edges bounding real cells )
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- if(iEdge <= nEdges) 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)
- enddo
- endif
+ end do
end do
end do
+
!
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
enddo
enddo
@@ -963,7 +1181,7 @@
!
do iEdge = 1,nEdges
do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
enddo
enddo
Deleted: branches/ocean_projects/topography2_mrp/src/core_sw/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_sw/mpas_interface.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/core_sw/mpas_interface.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,90 +0,0 @@
-subroutine mpas_setup_test_case(domain)
-
- use grid_types
- use test_cases
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- call setup_sw_test_case(domain)
-
-end subroutine mpas_setup_test_case
-
-
-subroutine mpas_init(block, mesh, dt)
-
- use grid_types
- use time_integration
- use RBF_interpolation
- use vector_reconstruction
-
- implicit none
-
- type (block_type), intent(inout) :: block
- type (mesh_type), intent(inout) :: mesh
- real (kind=RKIND), intent(in) :: dt
-
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(block % state % time_levs(1) % state, mesh)
-
-end subroutine mpas_init
-
-
-subroutine mpas_query(key, ivalue)
-
- implicit none
-
- character (len=256), intent(in) :: key
- integer, intent(out) :: ivalue
-
- if (index(key,'STORAGE_FACTOR') /= 0) then
- ivalue = 2
- end if
-
-end subroutine mpas_query
-
-
-subroutine mpas_timestep(domain, itimestep, dt)
-
- 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)
-
- if(config_stats_interval .gt. 0) then
- if(mod(itimestep, config_stats_interval) == 0) then
- 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 % state % time_levs(2) % state, block_ptr % mesh, &
- itimestep, dt)
- call timer_stop("global_diagnostics")
- end if
- end if
-
-end subroutine mpas_timestep
-
-
-subroutine mpas_finalize()
-
- implicit none
-
-end subroutine mpas_finalize
Modified: branches/ocean_projects/topography2_mrp/src/driver/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/src/driver/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/driver/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,13 +1,13 @@
.SUFFIXES: .F .o
-OBJS = module_subdriver.o \
+OBJS = module_mpas_subdriver.o \
mpas.o
all: $(OBJS)
-module_subdriver.o:
+module_mpas_subdriver.o:
-mpas.o: module_subdriver.o
+mpas.o: module_mpas_subdriver.o
clean:
$(RM) *.o *.mod *.f90
Copied: branches/ocean_projects/topography2_mrp/src/driver/module_mpas_subdriver.F (from rev 616, trunk/mpas/src/driver/module_mpas_subdriver.F)
===================================================================
--- branches/ocean_projects/topography2_mrp/src/driver/module_mpas_subdriver.F (rev 0)
+++ branches/ocean_projects/topography2_mrp/src/driver/module_mpas_subdriver.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -0,0 +1,86 @@
+module mpas_subdriver
+
+ use mpas_framework
+ use mpas_core
+
+ type (dm_info), pointer :: dminfo
+ type (domain_type), pointer :: domain
+ type (io_output_object) :: output_obj
+ integer :: output_frame
+
+
+ contains
+
+
+ subroutine mpas_init()
+
+ implicit none
+
+ real (kind=RKIND) :: dt
+
+ call timer_start("total time")
+ call timer_start("initialize")
+
+
+ !
+ ! Initialize infrastructure
+ !
+ call mpas_framework_init(dminfo, domain)
+
+
+ call input_state_for_domain(domain)
+
+
+ !
+ ! Initialize core
+ !
+ call mpas_core_init(domain)
+
+ call timer_stop("initialize")
+
+
+ !
+ ! Set up output streams to be written to by the MPAS core
+ !
+ output_frame = 1
+ call output_state_init(output_obj, domain, "OUTPUT")
+
+ end subroutine mpas_init
+
+
+ subroutine mpas_run()
+
+ implicit none
+
+ call mpas_core_run(domain, output_obj, output_frame)
+
+ end subroutine mpas_run
+
+
+ subroutine mpas_finalize()
+
+ implicit none
+
+ !
+ ! Finalize output streams
+ !
+ call output_state_finalize(output_obj, domain % dminfo)
+
+
+ !
+ ! Finalize core
+ !
+ call mpas_core_finalize(domain)
+
+ call timer_stop("total time")
+ call timer_write()
+
+
+ !
+ ! Finalize infrastructure
+ !
+ call mpas_framework_finalize(dminfo, domain)
+
+ end subroutine mpas_finalize
+
+end module mpas_subdriver
Deleted: branches/ocean_projects/topography2_mrp/src/driver/module_subdriver.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/driver/module_subdriver.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/driver/module_subdriver.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,150 +0,0 @@
-module subdriver
-
- use grid_types
- use io_output
- use configure
- use dmpar
- use timer
-
- integer :: output_frame
- integer :: restart_frame
-
- interface
- subroutine mpas_init(block, mesh, dt)
- use grid_types
- type (block_type), intent(inout) :: block
- type (mesh_type), intent(inout) :: mesh
- real (kind=RKIND), intent(in) :: dt
- end subroutine mpas_init
-
- subroutine mpas_timestep(domain, itimestep, dt)
- use grid_types
- type (domain_type), intent(inout) :: domain
- integer, intent(in) :: itimestep
- real (kind=RKIND), intent(in) :: dt
- end subroutine mpas_timestep
- end interface
-
-
- contains
-
-
- subroutine solve(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Solve 2d shallow water equations over the time range specified in the
- ! namelist, writing model state to an output file periodically
- !
- ! Input/Output: domain - grid metadata and model state, to be advanced forward
- ! in time
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: ntimesteps, itimestep
- real (kind=RKIND) :: dt
- type (block_type), pointer :: block_ptr
- type (io_output_object) :: output_obj
- type (io_output_object) :: restart_obj
-
-
- ! Eventually, dt should be domain specific
- dt = config_dt
- ntimesteps = config_ntimesteps
-
-
- ! Compute diagnostic fields needed in solve loop, and initialize
- ! simulation time to 0 for all blocks
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call mpas_init(block_ptr, block_ptr % mesh, dt)
- if (.not. config_do_restart) block_ptr % state % time_levs(1) % state % xtime % scalar = 0.0
- block_ptr => block_ptr % next
- end do
-
- ! Before integrating, write out the initial state
- output_frame = 1
- restart_frame = 1
- call output_state_init(output_obj, domain, "OUTPUT")
- call write_output_frame(output_obj, domain)
-
-
- ! During integration, time level 1 stores the model state at the beginning of the
- ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
- do itimestep = 1,ntimesteps
- write(0,*) 'Doing timestep ', itimestep
- call timer_start("time integration")
- call mpas_timestep(domain, itimestep, dt)
- call timer_stop("time integration")
-
- ! Move time level 2 fields back into time level 1 for next time step
- call shift_time_levels_state(domain % blocklist % state)
-
- if (mod(itimestep, config_output_interval) == 0) then
- call write_output_frame(output_obj, domain)
- end if
- if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
- if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
- call output_state_for_domain(restart_obj, domain, restart_frame)
- restart_frame = restart_frame + 1
- end if
- end do
-
- call output_state_finalize(output_obj, domain % dminfo)
- if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
-
- end subroutine solve
-
-
- subroutine write_output_frame(output_obj, domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields for a domain and write model state to output file
- !
- ! Input/Output: domain - contains model state; diagnostic field are computed
- ! before returning
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- type (io_output_object), intent(inout) :: output_obj
-
- integer :: i, j, k
- integer :: eoe
- type (block_type), pointer :: block_ptr
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
- block_ptr => block_ptr % next
- end do
-
- call output_state_for_domain(output_obj, domain, output_frame)
- output_frame = output_frame + 1
-
- end subroutine write_output_frame
-
-
- subroutine compute_output_diagnostics(state, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields for a domain
- !
- ! Input: state - contains model prognostic fields
- ! grid - contains grid metadata
- !
- ! Output: state - upon returning, diagnostic fields will have be computed
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: state
- type (mesh_type), intent(in) :: grid
-
- integer :: i, eoe
- integer :: iEdge, k
-
- end subroutine compute_output_diagnostics
-
-
-end module subdriver
Modified: branches/ocean_projects/topography2_mrp/src/driver/mpas.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/driver/mpas.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/driver/mpas.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -1,51 +1,15 @@
program mpas
- use grid_types
- use configure
- use io_input
- use dmpar
- use timer
- use subdriver
+ use mpas_subdriver
implicit none
- interface mpas_setup_test_case
- subroutine mpas_setup_test_case(domain)
- use grid_types
- type (domain_type), intent(inout) :: domain
- end subroutine mpas_setup_test_case
- end interface mpas_setup_test_case
+ call mpas_init()
+ call mpas_run()
- type (dm_info), pointer :: dminfo
- type (domain_type), pointer :: domain
+ call mpas_finalize()
- allocate(dminfo)
- call dmpar_init(dminfo)
-
- call read_namelist(dminfo)
-
- call timer_start("total time")
-
- call timer_start("initialize")
- call allocate_domain(domain, dminfo)
-
- call input_state_for_domain(domain)
-
- call mpas_init_domain(domain)
-
- if (.not. config_do_restart) call mpas_setup_test_case(domain)
- call timer_stop("initialize")
-
- call solve(domain)
-
- call deallocate_domain(domain)
-
- call timer_stop("total time")
- call timer_write()
-
- call dmpar_finalize(dminfo)
-
stop
end program mpas
Modified: branches/ocean_projects/topography2_mrp/src/framework/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/src/framework/Makefile 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/framework/Makefile 2010-11-15 23:08:12 UTC (rev 617)
@@ -4,7 +4,8 @@
ZOLTANOBJ = module_zoltan_interface.o
endif
-OBJS = module_timer.o \
+OBJS = module_mpas_framework.o \
+ module_timer.o \
module_configure.o \
module_constants.o \
module_grid_types.o \
@@ -22,6 +23,8 @@
framework: $(OBJS)
ar -ru libframework.a $(OBJS)
+module_mpas_framework.o: module_dmpar.o module_io_input.o module_io_output.o module_grid_types.o module_configure.o module_timer.o
+
module_configure.o: module_dmpar.o
module_grid_types.o: module_dmpar.o
Modified: branches/ocean_projects/topography2_mrp/src/framework/module_grid_types.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/framework/module_grid_types.F 2010-11-15 21:53:50 UTC (rev 616)
+++ branches/ocean_projects/topography2_mrp/src/framework/module_grid_types.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -93,7 +93,6 @@
! Derived type for storing part of a domain; used as a basic unit of work for a process
type block_type
- integer :: storageFactor ! Additional storage used by time integration scheme
#include "block_group_members.inc"
@@ -113,17 +112,7 @@
type (dm_info), pointer :: dminfo
end type domain_type
- !
- ! Solver interface routines provided by specific dycore
- !
- interface
- subroutine mpas_query(key, ivalue)
- character (len=256), intent(in) :: key
- integer, intent(out) :: ivalue
- end subroutine mpas_query
- end interface
-
contains
@@ -152,14 +141,10 @@
#include "dim_dummy_decls.inc"
integer :: i
- character (len=256) :: key
nullify(b % prev)
nullify(b % next)
- key = 'STORAGE_FACTOR'
- call mpas_query(key, b % storageFactor)
-
allocate(b % parinfo)
b % domain => dom
Copied: branches/ocean_projects/topography2_mrp/src/framework/module_mpas_framework.F (from rev 616, trunk/mpas/src/framework/module_mpas_framework.F)
===================================================================
--- branches/ocean_projects/topography2_mrp/src/framework/module_mpas_framework.F (rev 0)
+++ branches/ocean_projects/topography2_mrp/src/framework/module_mpas_framework.F 2010-11-15 23:08:12 UTC (rev 617)
@@ -0,0 +1,44 @@
+module mpas_framework
+
+ use dmpar
+ use grid_types
+ use io_input
+ use io_output
+ use configure
+ use timer
+
+
+ contains
+
+
+ subroutine mpas_framework_init(dminfo, domain)
+
+ implicit none
+
+ type (dm_info), pointer :: dminfo
+ type (domain_type), pointer :: domain
+
+ allocate(dminfo)
+ call dmpar_init(dminfo)
+
+ call read_namelist(dminfo)
+
+ call allocate_domain(domain, dminfo)
+
+ end subroutine mpas_framework_init
+
+
+ subroutine mpas_framework_finalize(dminfo, domain)
+
+ implicit none
+
+ type (dm_info), pointer :: dminfo
+ type (domain_type), pointer :: domain
+
+ call deallocate_domain(domain)
+
+ call dmpar_finalize(dminfo)
+
+ end subroutine mpas_framework_finalize
+
+end module mpas_framework
</font>
</pre>