<p><b>mpetersen@lanl.gov</b> 2012-05-22 10:27:19 -0600 (Tue, 22 May 2012)</p><p>BRANCH COMMIT: Computation of transport through sections within MPAS-Ocean is complete. The output file now contains the variable acc_uTransportSection ( nSections Time ), which contains a single value for each section, which is the transport across the section in Sv. There is a bit-for-bit match in the transport between this method and using the matlab code in mpas/branches/tools/transport_sections.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/transport_sections/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/transport_sections/src/core_ocean/Registry 2012-05-22 16:18:06 UTC (rev 1927)
+++ branches/ocean_projects/transport_sections/src/core_ocean/Registry 2012-05-22 16:27:19 UTC (rev 1928)
@@ -20,7 +20,7 @@
namelist integer io config_frames_per_outfile 0
namelist integer io config_pio_num_iotasks 0
namelist integer io config_pio_stride 1
-namelist integer io config_record_transport false
+namelist logical io config_record_transport false
namelist character decomposition config_block_decomp_file_prefix graph.info.part.
namelist integer decomposition config_number_of_blocks 0
namelist logical decomposition config_explicit_proc_decomp .false.
@@ -312,7 +312,7 @@
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
-var persistent real acc_uTransportSection ( maxNEdgesInSection nSections Time ) 2 - acc_uTransportSection state - -
+var persistent real acc_uTransportSection ( nSections Time ) 2 o acc_uTransportSection state - -
% Variables that contain sections for recording transport, read from grid or restart file:
var persistent integer sectionEdgeIndex ( maxNEdgesInSection nSections ) 0 - sectionEdgeIndex mesh - -
Modified: branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_mpas_core.F 2012-05-22 16:18:06 UTC (rev 1927)
+++ branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_mpas_core.F 2012-05-22 16:27:19 UTC (rev 1928)
@@ -232,6 +232,137 @@
end subroutine ocn_simulation_clock_init!}}}
+ subroutine ocn_initialize_transport_sections(mesh,err)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Read in transport section variables from text files.
+ ! This is temporary, will only be needed while these variables are not read
+ ! in properly from the netcdf file.
+ !
+ ! Input: mesh - contains mesh metadata
+ !
+ ! Output: mesh - upon return, section variables are filled.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use mpas_grid_types
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, intent(out) :: err
+ integer :: i, j
+ integer :: iEdge, k
+
+ err = 0
+
+ open (unit=15, file = "sectionEdgeIndex.txt")
+ read(15,*) mesh % sectionEdgeIndex % array
+ close(15)
+
+ open (unit=15, file = "sectionEdgeSign.txt")
+ read(15,*) mesh % sectionEdgeSign % array
+ close(15)
+
+ open (unit=15, file = "sectionCoord.txt")
+ read(15,*) mesh % sectionCoord % array
+ close(15)
+
+! mrp debugging output
+ print *, 'mesh % nSections',mesh % nSections
+ print *, 'mesh % maxNEdgesInSection',mesh % maxNEdgesInSection
+ print *, 'mesh % latLonPairs',mesh % latLonPairs
+ print *, 'mesh % CharLength8',mesh % CharLength8
+ print *, 'mesh % CharLength120',mesh % CharLength120
+
+ print *, 'size(mesh % nEdgesInSection',size(mesh % nEdgesInSection % array)
+ print *, 'size(mesh % sectionEdgeIndex',size(mesh % sectionEdgeIndex % array,1)
+ print *, 'size(mesh % sectionEdgeIndex',size(mesh % sectionEdgeIndex % array,2)
+ print *, 'size(mesh % sectionEdgeSign',size(mesh % sectionEdgeSign % array,1)
+ print *, 'size(mesh % sectionEdgeSign',size(mesh % sectionEdgeSign % array,2)
+ print *, 'size(mesh % sectionCoord',size(mesh % sectionCoord % array,1)
+ print *, 'size(mesh % sectionCoord',size(mesh % sectionCoord % array,2)
+ do j=1,mesh % nSections
+ print *, 'nSection = ',j
+ print *, 'mesh % nEdgesInSection',mesh % nEdgesInSection % array(j)
+ print *, 'mesh % sectionCoord',mesh % sectionCoord % array(:,j)
+ print *, 'mesh % sectionEdgeIndex',mesh % sectionEdgeIndex % array(:,j)
+ print *, 'mesh % sectionEdgeSign',mesh % sectionEdgeSign % array(:,j)
+ end do
+! mrp debugging output end
+
+ end subroutine ocn_initialize_transport_sections!}}}
+
+
+ subroutine ocn_find_local_transport_sections(mesh,err)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Read in transport section variables from text files.
+ ! This is temporary, will only be needed while these variables are not read
+ ! in properly from the netcdf file.
+ !
+ ! Input: mesh - contains mesh metadata
+ !
+ ! Output: mesh - upon return, local section variables are filled.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use mpas_grid_types
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, intent(out) :: err
+ integer :: iSection, iEdgeInSection,iEdge,iLocal, j
+
+ integer :: &
+ nSections, nEdgesSolve
+ integer, dimension(:), pointer :: &
+ nEdgesInSection, nLocalEdgesInSection, indexToEdgeID
+ integer, dimension(:,:), pointer :: &
+ sectionEdgeIndex, sectionEdgeSign, &
+ localSectionEdgeIndex, localSectionEdgeSign
+
+ nEdgesSolve = mesh % nEdgesSolve
+ nSections = mesh % nSections
+ nEdgesInSection => mesh % nEdgesInSection % array
+ nLocalEdgesInSection => mesh % nLocalEdgesInSection % array
+ indexToEdgeID => mesh % indexToEdgeID % array
+ sectionEdgeIndex => mesh % SectionEdgeIndex % array
+ localSectionEdgeIndex => mesh % LocalSectionEdgeIndex % array
+ sectionEdgeSign => mesh % SectionEdgeSign % array
+ localSectionEdgeSign => mesh % LocalSectionEdgeSign % array
+
+ err = 0
+
+ localSectionEdgeIndex = 0
+ localSectionEdgeSign = 0
+ do iSection = 1, nSections
+ iLocal = 0
+ do iEdgeInSection = 1,nEdgesInSection(iSection)
+ do iEdge = 1, nEdgesSolve
+ if (sectionEdgeIndex(iEdgeInSection,iSection).eq.indexToEdgeID(iEdge)) then
+ iLocal = iLocal + 1
+ localSectionEdgeIndex(iLocal,iSection) = iEdge
+ localSectionEdgeSign(iLocal,iSection) = sectionEdgeSign(iEdgeInSection,iSection)
+ exit
+ end if
+ end do
+ end do
+ nLocalEdgesInSection(iSection) = iLocal
+ end do
+
+! mrp debugging output
+ do j=1,mesh % nSections
+ print *, 'nSection = ',j
+ print *, 'mesh % nLocalEdgesInSection',mesh % nLocalEdgesInSection % array(j)
+ print *, 'mesh % localSectionEdgeIndex',mesh % localSectionEdgeIndex % array(:,j)
+ print *, 'mesh % localSectionEdgeSign',mesh % localSectionEdgeSign % array(:,j)
+ end do
+! mrp debugging output end
+
+ end subroutine ocn_find_local_transport_sections!}}}
+
+
+
subroutine mpas_init_block(block, mesh, dt, err)!{{{
use mpas_grid_types
@@ -249,32 +380,15 @@
integer :: i, iEdge, iCell, k
integer :: err1
-! mrp debug
-
- print *, 'block % mesh % nSections',block % mesh % nSections
- print *, 'block % mesh % maxNEdgesInSection',block % mesh % maxNEdgesInSection
- print *, 'block % mesh % latLonPairs',block % mesh % latLonPairs
- print *, 'block % mesh % CharLength8',block % mesh % CharLength8
- print *, 'block % mesh % CharLength120',block % mesh % CharLength120
-
- print *, 'size(block % mesh % nEdgesInSection',size(block % mesh % nEdgesInSection % array)
- print *, 'block % mesh % nEdgesInSection',block % mesh % nEdgesInSection % array
- print *, 'size(block % mesh % sectionEdgeIndex',size(block % mesh % sectionEdgeIndex % array,1)
- print *, 'size(block % mesh % sectionEdgeIndex',size(block % mesh % sectionEdgeIndex % array,2)
- print *, 'block % mesh % sectionEdgeIndex',block % mesh % sectionEdgeIndex % array
- print *, 'size(block % mesh % sectionEdgeSign',size(block % mesh % sectionEdgeSign % array,1)
- print *, 'size(block % mesh % sectionEdgeSign',size(block % mesh % sectionEdgeSign % array,2)
- print *, 'block % mesh % sectionEdgeSign',block % mesh % sectionEdgeSign % array
- print *, 'size(block % mesh % sectionCoord',size(block % mesh % sectionCoord % array,1)
- print *, 'size(block % mesh % sectionCoord',size(block % mesh % sectionCoord % array,2)
- print *, 'block % mesh % sectionCoord',block % mesh % sectionCoord % array
-
-! mrp debug end
-
call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh, err1)
err = ior(err, err1)
+ if (config_record_transport) then
+ call ocn_initialize_transport_sections(mesh, err)
+ call ocn_find_local_transport_sections(mesh, err)
+ end if
+
call ocn_time_average_init(block % state % time_levs(1) % state)
call mpas_timer_start("diagnostic solve", .false., initDiagSolveTimer)
@@ -416,7 +530,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call ocn_time_average_normalize(block_ptr % state % time_levs(1) % state)
+ call ocn_time_average_normalize(block_ptr % mesh,block_ptr % state % time_levs(1) % state,domain % dminfo)
block_ptr => block_ptr % next
end do
Modified: branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_average.F
===================================================================
--- branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_average.F 2012-05-22 16:18:06 UTC (rev 1927)
+++ branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_average.F 2012-05-22 16:27:19 UTC (rev 1928)
@@ -1,6 +1,8 @@
module ocn_time_average
+ use mpas_configure
use mpas_grid_types
+ use mpas_dmpar
implicit none
save
@@ -13,8 +15,8 @@
real, pointer :: nAccumulate
- real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar
- real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar, acc_uTransportSection
+ real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar, acc_uTransportSection
+ real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar
real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar
nAccumulate => state % nAccumulate % scalar
@@ -43,7 +45,8 @@
end subroutine ocn_time_average_init!}}}
- subroutine ocn_time_average_accumulate(state, old_state)!{{{
+ subroutine ocn_time_average_accumulate(mesh,state, old_state)!{{{
+ type (mesh_type), intent(in) :: mesh
type (state_type), intent(inout) :: state
type (state_type), intent(in) :: old_state
@@ -76,7 +79,6 @@
acc_uReconstructMeridionalVar => state % acc_uReconstructMeridionalVar % array
acc_u => state % acc_u % array
acc_uVar => state % acc_uVar % array
- acc_uTransportSection => state % acc_uTransportSection % array
old_acc_ssh => old_state % acc_ssh % array
old_acc_sshVar => old_state % acc_sshVar % array
@@ -86,7 +88,6 @@
old_acc_uReconstructMeridionalVar => old_state % acc_uReconstructMeridionalVar % array
old_acc_u => old_state % acc_u % array
old_acc_uVar => old_state % acc_uVar % array
- old_acc_uTransportSection => old_state % acc_uTransportSection % array
acc_ssh = old_acc_ssh + ssh
acc_sshVar = old_acc_sshVar + ssh**2
@@ -97,24 +98,78 @@
acc_u = old_acc_u + u
acc_uVar = old_acc_uVar + u**2
-! do iSection = 1,nSections
-! do i = 1,nLocalEdgesInSection(iSection)
-! acc_uTransportSection(i,iSection) &
-! = old_acc_uTransportSection(i,iSection) &
-! + u()
-! end do
-! end do
+ call ocn_transport_section_accumulate(mesh, state, old_state)
nAccumulate = old_nAccumulate + 1
end subroutine ocn_time_average_accumulate!}}}
- subroutine ocn_time_average_normalize(state)!{{{
+ subroutine ocn_transport_section_accumulate(mesh, state, old_state)!{{{
+ type (mesh_type), intent(in) :: mesh
type (state_type), intent(inout) :: state
+ type (state_type), intent(in) :: old_state
+ integer :: iSection, iLocalEdge, iEdge, k
+
+ integer :: &
+ nSections
+ integer, dimension(:), pointer :: &
+ nLocalEdgesInSection, maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: &
+ localSectionEdgeIndex, localSectionEdgeSign
+
+ real (kind=RKIND) :: tempSum, m3ps_to_Sv
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ acc_uTransportSection, old_acc_uTransportSection, dvEdge, hZLevel
+
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ u, h_edge
+
+ m3ps_to_Sv = 1e-6 ! m^3/sec flux to Sverdrups
+
+ nSections = mesh % nSections
+ nLocalEdgesInSection => mesh % nLocalEdgesInSection % array
+ maxLevelEdgeTop => mesh % maxLevelEdgeTop % array
+ dvEdge => mesh % dvEdge % array
+ localSectionEdgeIndex => mesh % localSectionEdgeIndex % array
+ localSectionEdgeSign => mesh % localSectionEdgeSign % array
+ hZLevel => mesh % hZLevel % array
+
+ u => state % u % array
+ h_edge => state % h_edge % array
+
+ acc_uTransportSection => state % acc_uTransportSection % array
+ old_acc_uTransportSection => old_state % acc_uTransportSection % array
+ do iSection = 1,nSections
+ tempSum = 0.0
+ do iLocalEdge = 1,nLocalEdgesInSection(iSection)
+ iEdge = localSectionEdgeIndex(iLocalEdge,iSection)
+ do k=1,maxLevelEdgeTop(iEdge)
+ tempSum = tempSum + &
+ + localSectionEdgeSign(iLocalEdge,iSection) &
+ * u(k,iEdge)*dvEdge(iEdge) &
+ * h_edge(k,iEdge)*m3ps_to_Sv;
+ end do
+ end do
+ acc_uTransportSection(iSection) &
+ = old_acc_uTransportSection(iSection) + tempSum
+ end do
+
+ end subroutine ocn_transport_section_accumulate!}}}
+
+ subroutine ocn_time_average_normalize(mesh,state,dminfo)!{{{
+ type (mesh_type), intent(in) :: mesh
+ type (state_type), intent(inout) :: state
+ type (dm_info) :: dminfo
+
+ integer :: &
+ nSections
+
real, pointer :: nAccumulate
- real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar
- real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar, acc_uTransportSection
+ real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar, acc_uTransportSection
+ real (kind=RKIND), dimension(:), allocatable :: temp_acc_uTransportSection
+ real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar
real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar
nAccumulate => state % nAccumulate % scalar
@@ -127,7 +182,6 @@
acc_uReconstructMeridionalVar => state % acc_uReconstructMeridionalVar % array
acc_u => state % acc_u % array
acc_uVar => state % acc_uVar % array
- acc_uTransportSection => state % acc_uTransportSection % array
acc_ssh = acc_ssh / nAccumulate
acc_sshVar = acc_sshVar / nAccumulate
@@ -137,8 +191,18 @@
acc_uReconstructMeridionalVar = acc_uReconstructMeridionalVar / nAccumulate
acc_u = acc_u / nAccumulate
acc_uVar = acc_uVar / nAccumulate
- acc_uTransportSection = acc_uTransportSection / nAccumulate
+ if (config_record_transport) then
+ nSections = mesh % nSections
+ acc_uTransportSection => state % acc_uTransportSection % array
+ allocate(temp_acc_uTransportSection(nSections))
+ temp_acc_uTransportSection = acc_uTransportSection / nAccumulate
+ call mpas_dmpar_sum_real_array(dminfo, nSections, temp_acc_uTransportSection, acc_uTransportSection)
+ deallocate(temp_acc_uTransportSection)
+ else
+ acc_uTransportSection = 0.0
+ end if
+
end subroutine ocn_time_average_normalize!}}}
end module ocn_time_average
Modified: branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-05-22 16:18:06 UTC (rev 1927)
+++ branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-05-22 16:27:19 UTC (rev 1928)
@@ -356,7 +356,7 @@
)
!TDR
- call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
+ call ocn_time_average_accumulate(block % mesh,block % state % time_levs(2) % state, block % state % time_levs(1) % state)
block => block % next
end do
Modified: branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_split.F 2012-05-22 16:18:06 UTC (rev 1927)
+++ branches/ocean_projects/transport_sections/src/core_ocean/mpas_ocn_time_integration_split.F 2012-05-22 16:27:19 UTC (rev 1928)
@@ -888,7 +888,7 @@
!TDR
- call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
+ call ocn_time_average_accumulate(block % mesh, block % state % time_levs(2) % state, block % state % time_levs(1) % state)
block => block % next
</font>
</pre>