From d4be0422c708ebf37aac0988166b7d045570385b Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 13 Apr 2017 22:58:16 -0600 Subject: [PATCH 01/11] layer volume average fixes --- ...egistry_layer_volume_weighted_averages.xml | 11 +- .../mpas_ocn_layer_volume_weighted_averages.F | 262 ++++++++++++------ 2 files changed, 181 insertions(+), 92 deletions(-) diff --git a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml index e76c2eb071..6296cd775d 100644 --- a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml +++ b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml @@ -31,7 +31,10 @@ - + + diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index 1099e7081a..996b49ba98 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -87,6 +87,13 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain + type (mpas_pool_type), pointer :: layerVolumeWeightedAverageAMPool + type (mpas_pool_type), pointer :: meshPool + real (kind=RKIND), dimension(:,:,:), pointer :: layerVolumeMask + integer, pointer :: nOceanRegionsTmp, nCellsSolve, nCells, nVertLevels + integer, pointer :: nLayerVolWeightedAvgFields + real (kind=RKIND), dimension(:), pointer :: lonCell, latCell + integer, dimension(:), pointer :: maxLevelCell !----------------------------------------------------------------- ! @@ -101,9 +108,96 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ ! local variables ! !----------------------------------------------------------------- + call mpas_pool_get_subpool(domain % blocklist % structs, 'layerVolumeWeightedAverageAM', layerVolumeWeightedAverageAMPool) + + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'layerVolumeMask', layerVolumeMask) + print *, '88234' + call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nOceanRegionsTmp', nOceanRegionsTmp) + call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nLayerVolWeightedAvgFields', nLayerVolWeightedAvgFields) + call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nCells', nCells) + call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool) + + call mpas_pool_get_array(meshPool, 'lonCell', lonCell) + call mpas_pool_get_array(meshPool, 'latCell', latCell) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) err = 0 + ! compute mask + call compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRegionsTmp, & + lonCell, latCell, layerVolumeMask) + + contains + + subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRegionsTmp, & + lonCell, latCell, workMask) + ! this subroutines produces a 0/1 mask that is multiplied with workArray to + ! allow for min/max/avg to represent specific regions of the ocean domain + ! + ! NOTE: computes_mask is temporary. workMask should be intent(in) to this entire module ! + ! + integer, intent(in) :: nVertLevels, nCells, nCellsSolve, nOceanRegionsTmp + integer, intent(in), dimension(:) :: maxLevelCell + real(kind=RKIND), dimension(:), intent(in) :: lonCell, latCell + real(kind=RKIND), dimension(:,:,:), intent(out) :: workMask + integer :: iCell, iLevel, iRegion + real(kind=RKIND) :: dtr + + dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND + workMask(:,:,:) = 0.0_RKIND + write(stderrUnit,*), 'blah1123213123123' + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + if(iLevel.le.maxLevelCell(iCell)) workMask(iCell,:,iLevel) = 1.0_RKIND + enddo + enddo + + write(stderrUnit,*), 'blah1' + + ! Arctic + do iCell=1,nCellsSolve + if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell,1,:) = 0.0_RKIND + enddo + write(stderrUnit,*), 'blah2' + ! Equatorial + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND + if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND + enddo + write(stderrUnit,*), 'blah3' + ! Southern Ocean + do iCell=1,nCellsSolve + if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell,3,:) = 0.0_RKIND + enddo + ! Nino 3 + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + enddo + ! Nino 4 + write(stderrUnit,*), 'blah4' + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + enddo + ! Nino 3.4 + write(stderrUnit,*), 'blah5' + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + enddo + + end subroutine compute_mask + + end subroutine ocn_init_layer_volume_weighted_averages!}}} !*********************************************************************** @@ -188,10 +282,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ real (kind=RKIND), dimension(:), pointer :: areaCell, lonCell, latCell ! scratch space - type(field2DReal), pointer :: workArrayField - real (kind=RKIND), dimension(:,:), pointer :: workArray - type(field1DReal), pointer :: workMaskField, workMinField, workMaxField, workSumField - real (kind=RKIND), dimension(:), pointer :: workMask, workMin, workMax, workSum + type(field2DReal), pointer :: workArrayField, workMinField, workMaxField, workSumField + real (kind=RKIND), dimension(:,:), pointer :: workArray, workMin, workMax, workSum + real (kind=RKIND), dimension(:,:,:), pointer :: layerVolumeMask ! local variables integer :: iDataField, nDefinedDataFields @@ -241,6 +334,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ maxValueWithinOceanVolumeRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'avgValueWithinOceanVolumeRegion', & avgValueWithinOceanVolumeRegion) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'layerVolumeMask', layerVolumeMask) ! loop over blocks ! NOTE: code is not valid for multiple blocks ! @@ -250,17 +344,17 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! get pointers to scratch variables call mpas_pool_get_subpool(block % structs, 'layerVolumeWeightedAverageAMScratch', scratchPool) call mpas_pool_get_field(scratchPool, 'workArrayLayerVolume', workArrayField) - call mpas_pool_get_field(scratchPool, 'workMaskLayerVolume', workMaskField) +! call mpas_pool_get_field(scratchPool, 'workMaskLayerVolume', workMaskField) call mpas_pool_get_field(scratchPool, 'workMinLayerVolume', workMinField) call mpas_pool_get_field(scratchPool, 'workMaxLayerVolume', workMaxField) call mpas_pool_get_field(scratchPool, 'workSumLayerVolume', workSumField) call mpas_allocate_scratch_field(workArrayField, .true.) - call mpas_allocate_scratch_field(workMaskField, .true.) +! call mpas_allocate_scratch_field(workMaskField, .true.) call mpas_allocate_scratch_field(workMinField, .true.) call mpas_allocate_scratch_field(workMaxField, .true.) call mpas_allocate_scratch_field(workSumField, .true.) workArray => workArrayField % array - workMask => workMaskField % array +! workMask => workMaskField % array workMin => workMinField % array workMax => workMaxField % array workSum => workSumField % array @@ -312,18 +406,15 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ workBufferMin(:) = +1.0e20_RKIND workBufferMax(:) = -1.0e20_RKIND - ! loop over all ocean regions - do iRegion=1,nOceanRegionsTmp - ! loop over the vertical do iLevel=1,nVertLevels ! compute mask - call compute_mask(iLevel, maxLevelCell, nCells, nCellsSolve, iRegion, lonCell, latCell, workMask) +! call compute_mask(iLevel, maxLevelCell, nCells, nCellsSolve, iRegion, lonCell, latCell, workMask) ! copy data into work array workArray( :,:) = 0.0_RKIND - workArray( 1,:) = workMask(:) +! workArray( 1,:) = workMask(:) workArray( 2,:) = areaCell(:) workArray( 3,:) = layerThickness(iLevel,:) workArray( 4,:) = density(iLevel,:) @@ -339,24 +430,26 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ workArray(14,:) = divergence(iLevel,:) workArray(15,:) = relativeVorticityCell(iLevel,:)*relativeVorticityCell(iLevel,:) - call compute_statistics(nDefinedDataFields, nCellsSolve, workArray, workMask, workMin, workMax, workSum) + call compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workArray, & + layerVolumeMask(:,:,iLevel), workMin, workMax, workSum) ! store data in buffer in order to allow only three dmpar calls - do iDataField=1,nDefinedDataFields - kBuffer = kBuffer+1 - workBufferSum(kBuffer) = workBufferSum(kBuffer) + workSum(iDataField) - workBufferMin(kBuffer) = min(workBufferMin(kBuffer), workMin(iDataField)) - workBufferMax(kBuffer) = max(workBufferMax(kBuffer), workMax(iDataField)) + do iRegion=1,nOceanRegionsTmp + do iDataField=1,nDefinedDataFields + kBuffer = kBuffer+1 + workBufferSum(kBuffer) = workBufferSum(kBuffer) + workSum(iRegion,iDataField) + workBufferMin(kBuffer) = min(workBufferMin(kBuffer), workMin(iRegion,iDataField)) + workBufferMax(kBuffer) = max(workBufferMax(kBuffer), workMax(iRegion,iDataField)) + enddo enddo enddo ! iLevel - end do ! iRegion kBuffer = 0 ! deallocate scratch fields call mpas_deallocate_scratch_field(workArrayField, .true.) - call mpas_deallocate_scratch_field(workMaskField, .true.) +! call mpas_deallocate_scratch_field(workMaskField, .true.) call mpas_deallocate_scratch_field(workMinField, .true.) call mpas_deallocate_scratch_field(workMaxField, .true.) call mpas_deallocate_scratch_field(workSumField, .true.) @@ -437,101 +530,94 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ contains - subroutine compute_mask(iLevel, maxLevelCell, nCells, nCellsSolve, iRegion, lonCell, latCell, workMask) + subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRegionsTmp, & + lonCell, latCell, workMask) ! this subroutines produces a 0/1 mask that is multiplied with workArray to ! allow for min/max/avg to represent specific regions of the ocean domain ! ! NOTE: computes_mask is temporary. workMask should be intent(in) to this entire module ! ! - integer, intent(in) :: iLevel, nCells, nCellsSolve, iRegion + integer, intent(in) :: nVertLevels, nCells, nCellsSolve, nOceanRegionsTmp integer, intent(in), dimension(:) :: maxLevelCell real(kind=RKIND), dimension(:), intent(in) :: lonCell, latCell - real(kind=RKIND), dimension(:), intent(out) :: workMask - integer :: iCell + real(kind=RKIND), dimension(:,:,:), intent(out) :: workMask + integer :: iCell, iLevel, iRegion real(kind=RKIND) :: dtr dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND - workMask(:) = 0.0_RKIND + workMask(:,:,:) = 0.0_RKIND do iCell=1,nCellsSolve - if(iLevel.le.maxLevelCell(iCell)) workMask(iCell) = 1.0_RKIND + do iLevel=1,nVertLevels + if(iLevel.le.maxLevelCell(iCell)) workMask(iCell,:,iLevel) = 1.0_RKIND + enddo enddo - if (iRegion.eq.1) then - ! Arctic - do iCell=1,nCellsSolve - if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - elseif (iRegion.eq.2) then - ! Equatorial - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - elseif (iRegion.eq.3) then - ! Southern Ocean - do iCell=1,nCellsSolve - if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - elseif (iRegion.eq.4) then - ! Nino 3 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - elseif (iRegion.eq.5) then - ! Nino 4 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - elseif (iRegion.eq.6) then - ! Nino 3.4 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(iCell) = 0.0_RKIND - enddo - else - ! global (do nothing!) - endif + ! Arctic + do iCell=1,nCellsSolve + if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell,1,:) = 0.0_RKIND + enddo + ! Equatorial + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND + if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND + enddo + ! Southern Ocean + do iCell=1,nCellsSolve + if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell,3,:) = 0.0_RKIND + enddo + ! Nino 3 + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + enddo + ! Nino 4 + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + enddo + ! Nino 3.4 + do iCell=1,nCellsSolve + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + enddo end subroutine compute_mask - subroutine compute_statistics(nDefinedDataFields, nCellsSolve, workArray, workMask, workMin, workMax, workSum) + subroutine compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workArray, workMask, workMin, workMax, workSum) ! this subroutines does the actual summing, min, max, masking ect ! this hides the messy code from the high-level subroutine - integer, intent(in) :: nDefinedDataFields, nCellsSolve - real(kind=RKIND), dimension(:,:), intent(in) :: workArray - real(kind=RKIND), dimension(:), intent(in) :: workMask - real(kind=RKIND), dimension(:), intent(out) :: workMin, workMax, workSum - integer :: iCell, iDataField + integer, intent(in) :: nDefinedDataFields, nCellsSolve, nOceanRegionsTmp + real(kind=RKIND), dimension(:,:), intent(in) :: workArray, workMask + real(kind=RKIND), dimension(:,:), intent(out) :: workMin, workMax, workSum + + integer :: iCell, iDataField,iRegion real(kind=RKIND) :: cellMask, cellArea, cellVolume workSum = 0.0_RKIND do iCell=1,nCellsSolve - cellMask = workMask(iCell) ! mask - cellArea = cellMask * workArray(2,iCell) ! area - cellVolume = cellArea * workArray(3,iCell) ! volume - workSum(1) = workSum(1) + cellMask ! 0/1 mask sum - workSum(2) = workSum(2) + cellArea ! area sum - workSum(3) = workSum(3) + cellVolume ! volume sum - do iDataField=4,nDefinedDataFields - workSum(iDataField) = workSum(iDataField) + cellVolume*workArray(iDataField,iCell) ! volume-weighted sum + do iRegion=1,nOceanRegionsTmp + cellMask = workMask(iCell, iRegion) ! mask + cellArea = cellMask * workArray(2,iCell) ! area + cellVolume = cellArea * workArray(3,iCell) ! volume + workSum(iRegion, 1) = workSum(iRegion, 1) + cellMask ! 0/1 mask sum + workSum(iRegion, 2) = workSum(iRegion, 2) + cellArea ! area sum + workSum(iRegion, 3) = workSum(iRegion, 3) + cellVolume ! volume sum + do iDataField=4,nDefinedDataFields + workSum(iRegion, iDataField) = workSum(iRegion, iDataField) + cellVolume*workArray(iDataField,iCell) ! volume-weighted sum + enddo + workMin(iRegion,iDatafield) = min(workSum(iRegion,iDataField), workMin(iRegion,iDataField)) + workMax(iRegion,iDataField) = max(workSum(iRegion,iDataField), workMax(iRegion,iDataField)) enddo enddo - do iDataField=1,nDefinedDataFields - workMin(iDataField) = minval(workArray(iDataField,1:nCellsSolve),workMask(1:nCellsSolve)>0.5_RKIND) - workMax(iDataField) = maxval(workArray(iDataField,1:nCellsSolve),workMask(1:nCellsSolve)>0.5_RKIND) - enddo - end subroutine compute_statistics end subroutine ocn_compute_layer_volume_weighted_averages!}}} From 568602843dcb837356607d288270eea2cebd20b3 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Fri, 14 Apr 2017 22:15:20 -0600 Subject: [PATCH 02/11] Cleans up unneeded code --- .../mpas_ocn_layer_volume_weighted_averages.F | 72 ------------------- 1 file changed, 72 deletions(-) diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index 996b49ba98..e06678dcb9 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -111,7 +111,6 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ call mpas_pool_get_subpool(domain % blocklist % structs, 'layerVolumeWeightedAverageAM', layerVolumeWeightedAverageAMPool) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'layerVolumeMask', layerVolumeMask) - print *, '88234' call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nOceanRegionsTmp', nOceanRegionsTmp) call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nLayerVolWeightedAvgFields', nLayerVolWeightedAvgFields) call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nVertLevels', nVertLevels) @@ -147,26 +146,21 @@ subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRe dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND workMask(:,:,:) = 0.0_RKIND - write(stderrUnit,*), 'blah1123213123123' do iCell=1,nCellsSolve do iLevel=1,nVertLevels if(iLevel.le.maxLevelCell(iCell)) workMask(iCell,:,iLevel) = 1.0_RKIND enddo enddo - write(stderrUnit,*), 'blah1' - ! Arctic do iCell=1,nCellsSolve if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell,1,:) = 0.0_RKIND enddo - write(stderrUnit,*), 'blah2' ! Equatorial do iCell=1,nCellsSolve if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND enddo - write(stderrUnit,*), 'blah3' ! Southern Ocean do iCell=1,nCellsSolve if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell,3,:) = 0.0_RKIND @@ -179,7 +173,6 @@ subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRe if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND enddo ! Nino 4 - write(stderrUnit,*), 'blah4' do iCell=1,nCellsSolve if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND @@ -187,7 +180,6 @@ subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRe if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND enddo ! Nino 3.4 - write(stderrUnit,*), 'blah5' do iCell=1,nCellsSolve if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND @@ -409,12 +401,8 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! loop over the vertical do iLevel=1,nVertLevels - ! compute mask -! call compute_mask(iLevel, maxLevelCell, nCells, nCellsSolve, iRegion, lonCell, latCell, workMask) - ! copy data into work array workArray( :,:) = 0.0_RKIND -! workArray( 1,:) = workMask(:) workArray( 2,:) = areaCell(:) workArray( 3,:) = layerThickness(iLevel,:) workArray( 4,:) = density(iLevel,:) @@ -530,66 +518,6 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ contains - subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRegionsTmp, & - lonCell, latCell, workMask) - ! this subroutines produces a 0/1 mask that is multiplied with workArray to - ! allow for min/max/avg to represent specific regions of the ocean domain - ! - ! NOTE: computes_mask is temporary. workMask should be intent(in) to this entire module ! - ! - integer, intent(in) :: nVertLevels, nCells, nCellsSolve, nOceanRegionsTmp - integer, intent(in), dimension(:) :: maxLevelCell - real(kind=RKIND), dimension(:), intent(in) :: lonCell, latCell - real(kind=RKIND), dimension(:,:,:), intent(out) :: workMask - integer :: iCell, iLevel, iRegion - real(kind=RKIND) :: dtr - - dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND - workMask(:,:,:) = 0.0_RKIND - do iCell=1,nCellsSolve - do iLevel=1,nVertLevels - if(iLevel.le.maxLevelCell(iCell)) workMask(iCell,:,iLevel) = 1.0_RKIND - enddo - enddo - - ! Arctic - do iCell=1,nCellsSolve - if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell,1,:) = 0.0_RKIND - enddo - ! Equatorial - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND - if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND - enddo - ! Southern Ocean - do iCell=1,nCellsSolve - if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell,3,:) = 0.0_RKIND - enddo - ! Nino 3 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - enddo - ! Nino 4 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - enddo - ! Nino 3.4 - do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - enddo - - end subroutine compute_mask - - subroutine compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workArray, workMask, workMin, workMax, workSum) ! this subroutines does the actual summing, min, max, masking ect ! this hides the messy code from the high-level subroutine From a184f5a0b858040e6998a1ef4af07ae6eb03fcb6 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 17 Apr 2017 10:15:07 -0600 Subject: [PATCH 03/11] Fixes array unpacking post send --- .../analysis_members/mpas_ocn_layer_volume_weighted_averages.F | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index e06678dcb9..5e7af829a9 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -437,7 +437,6 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! deallocate scratch fields call mpas_deallocate_scratch_field(workArrayField, .true.) -! call mpas_deallocate_scratch_field(workMaskField, .true.) call mpas_deallocate_scratch_field(workMinField, .true.) call mpas_deallocate_scratch_field(workMaxField, .true.) call mpas_deallocate_scratch_field(workSumField, .true.) @@ -452,8 +451,8 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! unpack the buffer into intent(out) of this analysis member kBuffer=0 + do iLevel=1,nVertLevels do iRegion=1,nOceanRegionsTmp - do iLevel=1,nVertLevels do iDataField=1,nDefinedDataFields kBuffer = kBuffer+1 avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion)=workBufferSumReduced(kBuffer) From 65839bc31343283b7c20a36fae473f728ba4c019 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 17 Apr 2017 13:59:12 -0600 Subject: [PATCH 04/11] Switch index ordering --- ...egistry_layer_volume_weighted_averages.xml | 21 ++-- .../mpas_ocn_layer_volume_weighted_averages.F | 111 +++++++++--------- 2 files changed, 61 insertions(+), 71 deletions(-) diff --git a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml index 6296cd775d..8dd1e78085 100644 --- a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml +++ b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml @@ -32,10 +32,10 @@ - - + @@ -129,7 +129,7 @@ description="Minimum relative enstrophy within region volume" /> - + @@ -223,7 +223,7 @@ description="Maximum relative enstrophy within region volume" /> - + @@ -319,13 +319,6 @@ - diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index 5e7af829a9..d03c54553a 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -87,7 +87,7 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain - type (mpas_pool_type), pointer :: layerVolumeWeightedAverageAMPool + type (mpas_pool_type), pointer :: layerVolumeWeightedAverageAMPool type (mpas_pool_type), pointer :: meshPool real (kind=RKIND), dimension(:,:,:), pointer :: layerVolumeMask integer, pointer :: nOceanRegionsTmp, nCellsSolve, nCells, nVertLevels @@ -146,45 +146,45 @@ subroutine compute_mask(nVertLevels, maxLevelCell, nCells, nCellsSolve, nOceanRe dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND workMask(:,:,:) = 0.0_RKIND - do iCell=1,nCellsSolve - do iLevel=1,nVertLevels - if(iLevel.le.maxLevelCell(iCell)) workMask(iCell,:,iLevel) = 1.0_RKIND + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + if(iLevel.le.maxLevelCell(iCell)) workMask(:,iCell,iLevel) = 1.0_RKIND enddo enddo ! Arctic do iCell=1,nCellsSolve - if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(iCell,1,:) = 0.0_RKIND + if(latCell(iCell).lt. 60.0_RKIND*dtr) workMask(1,iCell,:) = 0.0_RKIND enddo ! Equatorial do iCell=1,nCellsSolve - if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND - if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(iCell,2,:) = 0.0_RKIND + if(latCell(iCell).gt. 15.0_RKIND*dtr) workMask(2,iCell,:) = 0.0_RKIND + if(latCell(iCell).lt.-15.0_RKIND*dtr) workMask(2,iCell,:) = 0.0_RKIND enddo ! Southern Ocean do iCell=1,nCellsSolve - if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(iCell,3,:) = 0.0_RKIND + if(latCell(iCell).gt.-50.0_RKIND*dtr) workMask(3,iCell,:) = 0.0_RKIND enddo ! Nino 3 do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND - if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(iCell,4,:) = 0.0_RKIND + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(4,iCell,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(4,iCell,:) = 0.0_RKIND + if(lonCell(iCell).lt.210.0_RKIND*dtr) workMask(4,iCell,:) = 0.0_RKIND + if(lonCell(iCell).gt.270.0_RKIND*dtr) workMask(4,iCell,:) = 0.0_RKIND enddo ! Nino 4 do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND - if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(iCell,5,:) = 0.0_RKIND + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(5,iCell,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(5,iCell,:) = 0.0_RKIND + if(lonCell(iCell).lt.160.0_RKIND*dtr) workMask(5,iCell,:) = 0.0_RKIND + if(lonCell(iCell).gt.210.0_RKIND*dtr) workMask(5,iCell,:) = 0.0_RKIND enddo ! Nino 3.4 do iCell=1,nCellsSolve - if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND - if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(iCell,6,:) = 0.0_RKIND + if(latCell(iCell).gt. 5.0_RKIND*dtr) workMask(6,iCell,:) = 0.0_RKIND + if(latCell(iCell).lt. -5.0_RKIND*dtr) workMask(6,iCell,:) = 0.0_RKIND + if(lonCell(iCell).lt.190.0_RKIND*dtr) workMask(6,iCell,:) = 0.0_RKIND + if(lonCell(iCell).gt.240.0_RKIND*dtr) workMask(6,iCell,:) = 0.0_RKIND enddo end subroutine compute_mask @@ -336,17 +336,14 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! get pointers to scratch variables call mpas_pool_get_subpool(block % structs, 'layerVolumeWeightedAverageAMScratch', scratchPool) call mpas_pool_get_field(scratchPool, 'workArrayLayerVolume', workArrayField) -! call mpas_pool_get_field(scratchPool, 'workMaskLayerVolume', workMaskField) call mpas_pool_get_field(scratchPool, 'workMinLayerVolume', workMinField) call mpas_pool_get_field(scratchPool, 'workMaxLayerVolume', workMaxField) call mpas_pool_get_field(scratchPool, 'workSumLayerVolume', workSumField) call mpas_allocate_scratch_field(workArrayField, .true.) -! call mpas_allocate_scratch_field(workMaskField, .true.) call mpas_allocate_scratch_field(workMinField, .true.) call mpas_allocate_scratch_field(workMaxField, .true.) call mpas_allocate_scratch_field(workSumField, .true.) workArray => workArrayField % array -! workMask => workMaskField % array workMin => workMinField % array workMax => workMaxField % array workSum => workSumField % array @@ -425,9 +422,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iRegion=1,nOceanRegionsTmp do iDataField=1,nDefinedDataFields kBuffer = kBuffer+1 - workBufferSum(kBuffer) = workBufferSum(kBuffer) + workSum(iRegion,iDataField) - workBufferMin(kBuffer) = min(workBufferMin(kBuffer), workMin(iRegion,iDataField)) - workBufferMax(kBuffer) = max(workBufferMax(kBuffer), workMax(iRegion,iDataField)) + workBufferSum(kBuffer) = workBufferSum(kBuffer) + workSum(iDataField,iRegion) + workBufferMin(kBuffer) = min(workBufferMin(kBuffer), workMin(iDataField,iRegion)) + workBufferMax(kBuffer) = max(workBufferMax(kBuffer), workMax(iDataField,iRegion)) enddo enddo @@ -453,12 +450,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ kBuffer=0 do iLevel=1,nVertLevels do iRegion=1,nOceanRegionsTmp - do iDataField=1,nDefinedDataFields - kBuffer = kBuffer+1 - avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion)=workBufferSumReduced(kBuffer) - minValueWithinOceanLayerRegion(iDataField,iLevel,iRegion)=workBufferMinReduced(kBuffer) - maxValueWithinOceanLayerRegion(iDataField,iLevel,iRegion)=workBufferMaxReduced(kBuffer) - enddo + do iDataField=1,nDefinedDataFields + kBuffer = kBuffer+1 + avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferSumReduced(kBuffer) + minValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferMinReduced(kBuffer) + maxValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferMaxReduced(kBuffer) + enddo enddo enddo @@ -469,8 +466,8 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iRegion=1,nOceanRegionsTmp do iDataField=1,nDefinedDataFields do iLevel=1,nVertLevels - avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField, iRegion) & - + avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion) + avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField,iRegion) & + + avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) enddo enddo do iDataField=4,nDefinedDataFields @@ -488,25 +485,25 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! find min/max with region volume do iRegion=1,nOceanRegionsTmp do iDataField=1,nDefinedDataFields - minValueWithinOceanVolumeRegion(iDataField, iRegion) = minval(minValueWithinOceanLayerRegion(iDataField,:,iRegion)) - maxValueWithinOceanVolumeRegion(iDataField, iRegion) = maxval(minValueWithinOceanLayerRegion(iDataField,:,iRegion)) + minValueWithinOceanVolumeRegion(iDataField, iRegion) = minval(minValueWithinOceanLayerRegion(iDataField,iRegion,:)) + maxValueWithinOceanVolumeRegion(iDataField, iRegion) = maxval(minValueWithinOceanLayerRegion(iDataField,iRegion,:)) enddo enddo ! normalize averages layer-by-layer - do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - ! normalize all field by total volume in each layer - do iDataField=4,nDefinedDataFields - avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion) = avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion) & - / max(avgValueWithinOceanLayerRegion(3,iLevel,iRegion),1.0e-8_RKIND) - enddo - ! normalize total layer volume by layer area - avgValueWithinOceanLayerRegion(3,iLevel,iRegion) = avgValueWithinOceanLayerRegion(3,iLevel,iRegion) & - / max(avgValueWithinOceanLayerRegion(2,iLevel,iRegion),1.0e-8_RKIND) - ! normalize total layer area by number of cells in region - avgValueWithinOceanLayerRegion(2,iLevel,iRegion) = avgValueWithinOceanLayerRegion(2,iLevel,iRegion) & - / max(avgValueWithinOceanLayerRegion(1,iLevel,iRegion),1.0e-8_RKIND) + do iRegion=1,nOceanRegionsTmp + ! normalize all field by total volume in each layer + do iDataField=4,nDefinedDataFields + avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) = avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) & + / max(avgValueWithinOceanLayerRegion(3,iRegion,iLevel),1.0e-8_RKIND) + enddo + ! normalize total layer volume by layer area + avgValueWithinOceanLayerRegion(3,iRegion,iLevel) = avgValueWithinOceanLayerRegion(3,iRegion,iLevel) & + / max(avgValueWithinOceanLayerRegion(2,iRegion,iLevel),1.0e-8_RKIND) + ! normalize total layer area by number of cells in region + avgValueWithinOceanLayerRegion(2,iRegion,iLevel) = avgValueWithinOceanLayerRegion(2,iRegion,iLevel) & + / max(avgValueWithinOceanLayerRegion(1,iRegion,iLevel),1.0e-8_RKIND) enddo enddo @@ -531,18 +528,18 @@ subroutine compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workSum = 0.0_RKIND do iCell=1,nCellsSolve do iRegion=1,nOceanRegionsTmp - cellMask = workMask(iCell, iRegion) ! mask + cellMask = workMask(iRegion,iCell) ! mask cellArea = cellMask * workArray(2,iCell) ! area cellVolume = cellArea * workArray(3,iCell) ! volume - workSum(iRegion, 1) = workSum(iRegion, 1) + cellMask ! 0/1 mask sum - workSum(iRegion, 2) = workSum(iRegion, 2) + cellArea ! area sum - workSum(iRegion, 3) = workSum(iRegion, 3) + cellVolume ! volume sum - do iDataField=4,nDefinedDataFields - workSum(iRegion, iDataField) = workSum(iRegion, iDataField) + cellVolume*workArray(iDataField,iCell) ! volume-weighted sum + workSum(1,iRegion) = workSum(1,iRegion) + cellMask ! 0/1 mask sum + workSum(2,iRegion) = workSum(2,iRegion) + cellArea ! area sum + workSum(3,iRegion) = workSum(3,iRegion) + cellVolume ! volume sum + do iDataField=4,nDefinedDataFields + workSum(iDataField,iRegion) = workSum(iDataField,iRegion) + cellVolume*workArray(iDataField,iCell) ! volume-weighted sum + workMin(iDataField,iRegion) = min(workSum(iDataField,iRegion), workMin(iDataField,iRegion)) + workMax(iDataField,iRegion) = max(workSum(iDataField,iRegion), workMax(iDataField,iRegion)) + enddo enddo - workMin(iRegion,iDatafield) = min(workSum(iRegion,iDataField), workMin(iRegion,iDataField)) - workMax(iRegion,iDataField) = max(workSum(iRegion,iDataField), workMax(iRegion,iDataField)) - enddo enddo end subroutine compute_statistics From 645984b02c770c6b75edd196c3d12f1094b6c01c Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Wed, 19 Apr 2017 09:37:09 -0600 Subject: [PATCH 05/11] Add layer and min/max/avg to test suite comparison --- .../QU240/analysis_test/config_driver.xml | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/testing_and_setup/compass/ocean/global_ocean/QU240/analysis_test/config_driver.xml b/testing_and_setup/compass/ocean/global_ocean/QU240/analysis_test/config_driver.xml index 2ef14c3a98..d1cc01dcdf 100644 --- a/testing_and_setup/compass/ocean/global_ocean/QU240/analysis_test/config_driver.xml +++ b/testing_and_setup/compass/ocean/global_ocean/QU240/analysis_test/config_driver.xml @@ -36,8 +36,18 @@ + + - + + + + + + + + + From 1f66c6ceffc6bcab3f04ed10194799434cfbfd71 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Tue, 2 May 2017 15:48:04 -0600 Subject: [PATCH 06/11] Phil performance version 1: column ordering --- ...egistry_layer_volume_weighted_averages.xml | 45 +- .../mpas_ocn_layer_volume_weighted_averages.F | 934 ++++++++++++------ 2 files changed, 628 insertions(+), 351 deletions(-) diff --git a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml index 8dd1e78085..1f3aa5137b 100644 --- a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml +++ b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml @@ -1,8 +1,8 @@ - - @@ -31,11 +31,8 @@ - - - + + @@ -129,7 +126,7 @@ description="Minimum relative enstrophy within region volume" /> - + @@ -223,7 +220,7 @@ description="Maximum relative enstrophy within region volume" /> - + @@ -318,36 +315,6 @@ /> - - - - - - domain % blocklist + do while (associated(block)) + + numBlocks = numBlocks + 1 + call mpas_pool_get_dimension(block % dimensions, 'nCells', & + nCells) + maxCells = max(maxCells, nCells) + + block => block % next + end do + + ! allocate the region mask array + allocate(regionMask(maxCells,nRegions,numBlocks)) + + ! compute the region mask for each block + iBlock = 0 + block => domain % blocklist + do while (associated(block)) + iBlock = iBlock + 1 + + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', & + nCellsSolve) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'lonCell', lonCell) + call mpas_pool_get_array(meshPool, 'latCell', latCell) + + call ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & + lonCell, latCell) + + block => block % next + end do + + !----------------------------------------------------------------- end subroutine ocn_init_layer_volume_weighted_averages!}}} @@ -238,6 +222,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ !----------------------------------------------------------------- type (dm_info) :: dminfo + type (block_type), pointer :: block type (mpas_pool_type), pointer :: layerVolumeWeightedAverageAMPool type (mpas_pool_type), pointer :: statePool @@ -247,12 +232,11 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: tracersPool - real (kind=RKIND), dimension(:,:,:), pointer :: minValueWithinOceanLayerRegion - real (kind=RKIND), dimension(:,:,:), pointer :: maxValueWithinOceanLayerRegion - real (kind=RKIND), dimension(:,:,:), pointer :: avgValueWithinOceanLayerRegion - real (kind=RKIND), dimension(:,:), pointer :: minValueWithinOceanVolumeRegion - real (kind=RKIND), dimension(:,:), pointer :: maxValueWithinOceanVolumeRegion - real (kind=RKIND), dimension(:,:), pointer :: avgValueWithinOceanVolumeRegion + ! pointers to result arrays + real (kind=RKIND), dimension(:,:,:), pointer :: & + avgLayer, minLayer, maxLayer + real (kind=RKIND), dimension(:,:), pointer :: & + avgField, minField, maxField ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -268,25 +252,29 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ real (kind=RKIND), dimension(:,:), pointer :: divergence ! pointers to data in mesh pool - integer, pointer :: nVertLevels, nCells, nCellsSolve, nLayerVolWeightedAvgFields, nOceanRegionsTmp + integer, pointer :: nVertLevels, nCells, nCellsSolve, & + nFields, nRegions integer, pointer :: index_temperature, index_salinity integer, dimension(:), pointer :: maxLevelCell - real (kind=RKIND), dimension(:), pointer :: areaCell, lonCell, latCell - - ! scratch space - type(field2DReal), pointer :: workArrayField, workMinField, workMaxField, workSumField - real (kind=RKIND), dimension(:,:), pointer :: workArray, workMin, workMax, workSum - real (kind=RKIND), dimension(:,:,:), pointer :: layerVolumeMask + real (kind=RKIND), dimension(:), pointer :: areaCell ! local variables - integer :: iDataField, nDefinedDataFields - integer :: iCell, iLevel, iRegion, iTracer, err_tmp + integer :: iCell, iLevel, iRegion, iTracer, iField, iBlock + + real (kind=RKIND) :: maskTmp, normFactor ! for normalizing averages + + real (kind=RKIND), dimension(:,:), allocatable :: & + vertMask, &!mask for active vertical levels + fieldTmp !field array in local order ! buffers data for message passaging - integer :: kBuffer, kBufferLength - real (kind=RKIND), dimension(:), allocatable :: workBufferSum, workBufferSumReduced - real (kind=RKIND), dimension(:), allocatable :: workBufferMin, workBufferMinReduced - real (kind=RKIND), dimension(:), allocatable :: workBufferMax, workBufferMaxReduced + integer :: bufferAddr, bufferLength + real (kind=RKIND), dimension(:), allocatable :: & + msgBuffer, msgBufferReduced + + ! end of preamble + !----------------------------------------------------------------- + ! begin code ! assume no error err = 0 @@ -294,255 +282,357 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! set highest level pointer dminfo = domain % dminfo - ! find the number of regions, number of data fields and number of vertical levels - call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nOceanRegionsTmp', nOceanRegionsTmp) - call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nLayerVolWeightedAvgFields', nLayerVolWeightedAvgFields) - call mpas_pool_get_dimension(domain % blocklist % dimensions, 'nVertLevels', nVertLevels) - - ! allocate buffer for message passing - kBuffer=0 - kBufferLength=nOceanRegionsTmp*nLayerVolWeightedAvgFields*nVertLevels - allocate(workBufferSum(kBufferLength)) - allocate(workBufferMin(kBufferLength)) - allocate(workBufferMax(kBufferLength)) - allocate(workBufferSumReduced(kBufferLength)) - allocate(workBufferMinReduced(kBufferLength)) - allocate(workBufferMaxReduced(kBufferLength)) - workBufferSum=0.0_RKIND - workBufferMin=0.0_RKIND - workBufferMax=0.0_RKIND - workBufferSumReduced=0.0_RKIND - workBufferMinReduced=0.0_RKIND - workBufferMaxReduced=0.0_RKIND - - ! get pointers to analysis member arrays - call mpas_pool_get_subpool(domain % blocklist % structs, 'layerVolumeWeightedAverageAM', layerVolumeWeightedAverageAMPool) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'minValueWithinOceanLayerRegion', minValueWithinOceanLayerRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'maxValueWithinOceanLayerRegion', maxValueWithinOceanLayerRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'avgValueWithinOceanLayerRegion', avgValueWithinOceanLayerRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'minValueWithinOceanVolumeRegion', & - minValueWithinOceanVolumeRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'maxValueWithinOceanVolumeRegion', & - maxValueWithinOceanVolumeRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'avgValueWithinOceanVolumeRegion', & - avgValueWithinOceanVolumeRegion) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, 'layerVolumeMask', layerVolumeMask) + ! find the number of regions, fields and vertical levels + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nOceanRegionsTmp', nRegions) + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nLayerVolWeightedAvgFields', nFields) + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nVertLevels', nVertLevels) + + ! get pointers to analysis member arrays to store results + call mpas_pool_get_subpool(domain % blocklist % structs, & + 'layerVolumeWeightedAverageAM', & + layerVolumeWeightedAverageAMPool) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanLayerRegion', & + minLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanLayerRegion', & + maxLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanLayerRegion', & + avgLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanVolumeRegion', & + minField) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanVolumeRegion', & + maxField) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanVolumeRegion', & + avgField) + + ! initialize so this routine works for multiple blocks + avgLayer = 0.0_RKIND + minLayer = 1.e20_RKIND + maxLayer = -1.e20_RKIND ! loop over blocks - ! NOTE: code is not valid for multiple blocks ! + iBlock = 0 block => domain % blocklist do while (associated(block)) - ! get pointers to scratch variables - call mpas_pool_get_subpool(block % structs, 'layerVolumeWeightedAverageAMScratch', scratchPool) - call mpas_pool_get_field(scratchPool, 'workArrayLayerVolume', workArrayField) - call mpas_pool_get_field(scratchPool, 'workMinLayerVolume', workMinField) - call mpas_pool_get_field(scratchPool, 'workMaxLayerVolume', workMaxField) - call mpas_pool_get_field(scratchPool, 'workSumLayerVolume', workSumField) - call mpas_allocate_scratch_field(workArrayField, .true.) - call mpas_allocate_scratch_field(workMinField, .true.) - call mpas_allocate_scratch_field(workMaxField, .true.) - call mpas_allocate_scratch_field(workSumField, .true.) - workArray => workArrayField % array - workMin => workMinField % array - workMax => workMaxField % array - workSum => workSumField % array + iBlock = iBlock + 1 ! get pointers to pools - call mpas_pool_get_subpool(block % structs, 'state', statePool) - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) - call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool) - call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + call mpas_pool_get_subpool(block % structs, 'state', & + statePool) + call mpas_pool_get_subpool(block % structs, 'mesh', & + meshPool) + call mpas_pool_get_subpool(block % structs, 'diagnostics', & + diagnosticsPool) + call mpas_pool_get_subpool(block % structs, 'forcing', & + forcingPool) + call mpas_pool_get_subpool(statePool, 'tracers', & + tracersPool) ! get pointers to mesh - call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) - call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) - call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels) - call mpas_pool_get_dimension(block % dimensions, 'nLayerVolWeightedAvgFields', nLayerVolWeightedAvgFields) - call mpas_pool_get_dimension(block % dimensions, 'nOceanRegionsTmp', nOceanRegionsTmp) - call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature) - call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity) - call mpas_pool_get_array(meshPool, 'areaCell', areaCell) - call mpas_pool_get_array(meshPool, 'lonCell', lonCell) - call mpas_pool_get_array(meshPool, 'latCell', latCell) + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', & + nCellsSolve) + call mpas_pool_get_dimension(block % dimensions, 'nCells', & + nCells) + call mpas_pool_get_dimension(tracersPool, 'index_temperature', & + index_temperature) + call mpas_pool_get_dimension(tracersPool, 'index_salinity', & + index_salinity) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - ! test to make sure the arrays are big enough - nDefinedDataFields = size(avgValueWithinOceanLayerRegion,dim=1) - if (nDefinedDataFields > nLayerVolWeightedAvgFields) then - write (stderrUnit,*) 'Abort: nDefinedDataFields > nLayerVolWeightedAvgFields' - write (stderrUnit,*) ' : increase size of ocn_layer_volume_weighted_averages scratch space' - call mpas_dmpar_global_abort('MPAS-ocean: Abort: nDefinedDataFields > nLayerVolWeightedAvgFields') - endif - ! get pointers to data that will be analyzed - ! listed in the order in which the fields appear in {avg,min,max}ValueWithinOceanLayerRegion - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) - call mpas_pool_get_array(diagnosticsPool, 'density', density) - call mpas_pool_get_array(diagnosticsPool, 'potentialDensity', potentialDensity) - call mpas_pool_get_array(diagnosticsPool, 'BruntVaisalaFreqTop', BruntVaisalaFreqTop) - call mpas_pool_get_array(diagnosticsPool, 'velocityZonal', velocityZonal) - call mpas_pool_get_array(diagnosticsPool, 'velocityMeridional', velocityMeridional) - call mpas_pool_get_array(diagnosticsPool, 'vertVelocityTop', vertVelocityTop) - call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', kineticEnergyCell) - call mpas_pool_get_array(diagnosticsPool, 'relativeVorticityCell', relativeVorticityCell) - call mpas_pool_get_array(diagnosticsPool, 'divergence', divergence) - call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1) - - ! initialize buffers - workBufferSum(:) = 0.0_RKIND - workBufferMin(:) = +1.0e20_RKIND - workBufferMax(:) = -1.0e20_RKIND - - ! loop over the vertical + ! listed in the order in which the fields appear in output + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThickness, 1) + call mpas_pool_get_array(diagnosticsPool, 'density', & + density) + call mpas_pool_get_array(diagnosticsPool, 'potentialDensity', & + potentialDensity) + call mpas_pool_get_array(diagnosticsPool, 'BruntVaisalaFreqTop', & + BruntVaisalaFreqTop) + call mpas_pool_get_array(diagnosticsPool, 'velocityZonal', & + velocityZonal) + call mpas_pool_get_array(diagnosticsPool, 'velocityMeridional', & + velocityMeridional) + call mpas_pool_get_array(diagnosticsPool, 'vertVelocityTop', & + vertVelocityTop) + call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', & + kineticEnergyCell) + call mpas_pool_get_array(diagnosticsPool, 'relativeVorticityCell', & + relativeVorticityCell) + call mpas_pool_get_array(diagnosticsPool, 'divergence', & + divergence) + call mpas_pool_get_array(tracersPool, 'activeTracers', & + activeTracers, 1) + + + ! allocate arrays to hold local block masks, and fields + allocate(fieldTmp(nCellsSolve,nFields)) + !$omp single + allocate(vertMask(nCells,nVertLevels)) + !$omp end single + + call ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & + vertMask) + + !$omp do private(iLevel,iCell,iField,iRegion, & + !$omp maskTmp, fieldTmp) do iLevel=1,nVertLevels - ! copy data into work array - workArray( :,:) = 0.0_RKIND - workArray( 2,:) = areaCell(:) - workArray( 3,:) = layerThickness(iLevel,:) - workArray( 4,:) = density(iLevel,:) - workArray( 5,:) = potentialDensity(iLevel,:) - workArray( 6,:) = BruntVaisalaFreqTop(iLevel,:) - workArray( 7,:) = velocityZonal(iLevel,:) - workArray( 8,:) = velocityMeridional(iLevel,:) - workArray( 9,:) = vertVelocityTop(iLevel,:) - if ( associated(activeTracers) ) workArray(10,:) = activeTracers(index_temperature,iLevel,:) - if ( associated(activeTracers) ) workArray(11,:) = activeTracers(index_salinity,iLevel,:) - workArray(12,:) = kineticEnergyCell(iLevel,:) - workArray(13,:) = relativeVorticityCell(iLevel,:) - workArray(14,:) = divergence(iLevel,:) - workArray(15,:) = relativeVorticityCell(iLevel,:)*relativeVorticityCell(iLevel,:) - - call compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workArray, & - layerVolumeMask(:,:,iLevel), workMin, workMax, workSum) - - ! store data in buffer in order to allow only three dmpar calls - do iRegion=1,nOceanRegionsTmp - do iDataField=1,nDefinedDataFields - kBuffer = kBuffer+1 - workBufferSum(kBuffer) = workBufferSum(kBuffer) + workSum(iDataField,iRegion) - workBufferMin(kBuffer) = min(workBufferMin(kBuffer), workMin(iDataField,iRegion)) - workBufferMax(kBuffer) = max(workBufferMax(kBuffer), workMax(iDataField,iRegion)) - enddo - enddo - - enddo ! iLevel - - kBuffer = 0 - - ! deallocate scratch fields - call mpas_deallocate_scratch_field(workArrayField, .true.) - call mpas_deallocate_scratch_field(workMinField, .true.) - call mpas_deallocate_scratch_field(workMaxField, .true.) - call mpas_deallocate_scratch_field(workSumField, .true.) - + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask (iCell,iLevel) + fieldTmp(iCell,2) = areaCell (iCell) + fieldTmp(iCell,3) = layerThickness (iLevel,iCell) + fieldTmp(iCell,4) = density (iLevel,iCell) + fieldTmp(iCell,5) = potentialDensity (iLevel,iCell) + fieldTmp(iCell,6) = BruntVaisalaFreqTop (iLevel,iCell) + fieldTmp(iCell,7) = velocityZonal (iLevel,iCell) + fieldTmp(iCell,8) = velocityMeridional (iLevel,iCell) + fieldTmp(iCell,9) = vertVelocityTop (iLevel,iCell) + fieldTmp(iCell,10) = & + activeTracers(index_temperature,iLevel,iCell) + fieldTmp(iCell,11) = & + activeTracers(index_salinity,iLevel,iCell) + fieldTmp(iCell,12) = kineticEnergyCell (iLevel,iCell) + fieldTmp(iCell,13) = relativeVorticityCell(iLevel,iCell) + fieldTmp(iCell,14) = divergence (iLevel,iCell) + fieldTmp(iCell,15) = relativeVorticityCell(iLevel,iCell)* & + relativeVorticityCell(iLevel,iCell) + end do + + ! compute min and max + do iRegion=1,nRegions + do iField=1,nFields + do iCell=1,nCellsSolve + maskTmp = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + if (maskTmp /= 0.0_RKIND) then + minLayer(iField,iLevel,iRegion) = min( & + minLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + maxLayer(iField,iLevel,iRegion) = max( & + maxLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + endif + end do + end do + end do + + ! for the latter fields, we want volume-weighted sums + ! so replace thickness with cell volume + do iCell=1,nCellsSolve + fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) + end do + + ! compute layer averages for each region + do iRegion=1,nRegions + + ! replace field 1 with combined vertical, region mask + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + avgLayer(1,iLevel,iRegion) = & + avgLayer(1,iLevel,iRegion) + fieldTmp(iCell,1) + end do + + ! cell area + do iCell=1,nCellsSolve + avgLayer(2,iLevel,iRegion) = & + avgLayer(2,iLevel,iRegion) + fieldTmp(iCell,2)* & + fieldTmp(iCell,1) + end do + + ! for field three, we want cell volume instead of thickness + do iCell=1,nCellsSolve + avgLayer(3,iLevel,iRegion) = & + avgLayer(3,iLevel,iRegion) + fieldTmp(iCell,3)*& + fieldTmp(iCell,1) + end do + + ! for all remaining fields, we want a volume weighted avg + do iField=4,nFields + do iCell=1,nCellsSolve + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) + & + fieldTmp(iCell,iField)* & + fieldTmp(iCell,3) * & + fieldTmp(iCell,1) + end do + end do + end do ! iRegion + + end do !iLevel + + ! deallocate fields from this block + !$omp single + deallocate(fieldTmp) + deallocate(vertMask) + !$omp end single block => block % next end do - ! communication - call mpas_dmpar_sum_real_array(dminfo, kBufferLength, workBufferSum, workBufferSumReduced ) - call mpas_dmpar_min_real_array(dminfo, kBufferLength, workBufferMin, workBufferMinReduced ) - call mpas_dmpar_max_real_array(dminfo, kBufferLength, workBufferMax, workBufferMaxReduced ) + ! allocate buffer for message passing + bufferAddr=0 + bufferLength=nRegions*nFields*nVertLevels + !$omp single + allocate(msgBuffer (bufferLength)) + allocate(msgBufferReduced(bufferLength)) + !$omp end single + + ! communication for sums + ! pack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = avgLayer(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_sum_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + + ! compute normalized averages - ! unpack the buffer into intent(out) of this analysis member - kBuffer=0 + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions do iLevel=1,nVertLevels - do iRegion=1,nOceanRegionsTmp - do iDataField=1,nDefinedDataFields - kBuffer = kBuffer+1 - avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferSumReduced(kBuffer) - minValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferMinReduced(kBuffer) - maxValueWithinOceanLayerRegion(iDataField,iRegion,iLevel)=workBufferMaxReduced(kBuffer) - enddo - enddo + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + avgLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo enddo + enddo + !$omp end do ! compute vertical sum before layer-by-layer normalization - minValueWithinOceanVolumeRegion = 0.0_RKIND - maxValueWithinOceanVolumeRegion = 0.0_RKIND - avgValueWithinOceanVolumeRegion = 0.0_RKIND - do iRegion=1,nOceanRegionsTmp - do iDataField=1,nDefinedDataFields + !$omp do private(iLevel,iField,iRegion,normFactor) + do iRegion=1,nRegions + do iField=1,nFields + avgField(iField,iRegion) = 0.0_RKIND do iLevel=1,nVertLevels - avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField,iRegion) & - + avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) + avgField(iField,iRegion) = avgField(iField,iRegion) & + + avgLayer(iField,iLevel,iRegion) enddo enddo - do iDataField=4,nDefinedDataFields - avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField, iRegion) & - / max(avgValueWithinOceanVolumeRegion(3,iRegion),1.0e-8_RKIND) + normFactor = 1.0_RKIND / max(avgField(3,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgField(iField,iRegion) = avgField(iField,iRegion)*normFactor enddo ! normalize total region volume by total volume cell area - avgValueWithinOceanVolumeRegion(3,iRegion) = avgValueWithinOceanVolumeRegion(3,iRegion) & - / max(avgValueWithinOceanVolumeRegion(2,iRegion),1.0e-8_RKIND) + avgField(3,iRegion) = avgField(3,iRegion) & + / max(avgField(2,iRegion),1.0e-8_RKIND) ! normalize total volume cell area by total number of cells - avgValueWithinOceanVolumeRegion(2,iRegion) = avgValueWithinOceanVolumeRegion(2,iRegion) & - / max(avgValueWithinOceanVolumeRegion(1,iRegion),1.0e-8_RKIND) + avgField(2,iRegion) = avgField(2,iRegion) & + / max(avgField(1,iRegion),1.0e-8_RKIND) enddo + !$omp end do - ! find min/max with region volume - do iRegion=1,nOceanRegionsTmp - do iDataField=1,nDefinedDataFields - minValueWithinOceanVolumeRegion(iDataField, iRegion) = minval(minValueWithinOceanLayerRegion(iDataField,iRegion,:)) - maxValueWithinOceanVolumeRegion(iDataField, iRegion) = maxval(minValueWithinOceanLayerRegion(iDataField,iRegion,:)) - enddo + ! normalize averages layer-by-layer + !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + ! normalize all field by total volume in each layer + normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) * normFactor + enddo + ! normalize total layer volume by layer area + avgLayer(3,iLevel,iRegion) = avgLayer(3,iLevel,iRegion) & + / max(avgLayer(2,iLevel,iRegion),1.0e-8_RKIND) + ! normalize total layer area by number of cells in region + avgLayer(2,iLevel,iRegion) = avgLayer(2,iLevel,iRegion) & + / max(avgLayer(1,iLevel,iRegion),1.0e-8_RKIND) + enddo enddo + !$omp end do - ! normalize averages layer-by-layer + ! now compute min/max for each field, region + ! find global min/max + + ! pack the buffer for minval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions do iLevel=1,nVertLevels - do iRegion=1,nOceanRegionsTmp - ! normalize all field by total volume in each layer - do iDataField=4,nDefinedDataFields - avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) = avgValueWithinOceanLayerRegion(iDataField,iRegion,iLevel) & - / max(avgValueWithinOceanLayerRegion(3,iRegion,iLevel),1.0e-8_RKIND) - enddo - ! normalize total layer volume by layer area - avgValueWithinOceanLayerRegion(3,iRegion,iLevel) = avgValueWithinOceanLayerRegion(3,iRegion,iLevel) & - / max(avgValueWithinOceanLayerRegion(2,iRegion,iLevel),1.0e-8_RKIND) - ! normalize total layer area by number of cells in region - avgValueWithinOceanLayerRegion(2,iRegion,iLevel) = avgValueWithinOceanLayerRegion(2,iRegion,iLevel) & - / max(avgValueWithinOceanLayerRegion(1,iRegion,iLevel),1.0e-8_RKIND) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = minLayer(iField,iLevel,iRegion) + enddo enddo enddo + !$omp end do + ! communicate for global min + call mpas_dmpar_min_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + minLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do - ! deallocate buffers - deallocate(workBufferSumReduced) - deallocate(workBufferMinReduced) - deallocate(workBufferMaxReduced) - - contains - - subroutine compute_statistics(nDefinedDataFields, nOceanRegionsTmp, nCellsSolve, workArray, workMask, workMin, workMax, workSum) - ! this subroutines does the actual summing, min, max, masking ect - ! this hides the messy code from the high-level subroutine - - integer, intent(in) :: nDefinedDataFields, nCellsSolve, nOceanRegionsTmp - real(kind=RKIND), dimension(:,:), intent(in) :: workArray, workMask - real(kind=RKIND), dimension(:,:), intent(out) :: workMin, workMax, workSum - - integer :: iCell, iDataField,iRegion - real(kind=RKIND) :: cellMask, cellArea, cellVolume - - workSum = 0.0_RKIND - do iCell=1,nCellsSolve - do iRegion=1,nOceanRegionsTmp - cellMask = workMask(iRegion,iCell) ! mask - cellArea = cellMask * workArray(2,iCell) ! area - cellVolume = cellArea * workArray(3,iCell) ! volume - workSum(1,iRegion) = workSum(1,iRegion) + cellMask ! 0/1 mask sum - workSum(2,iRegion) = workSum(2,iRegion) + cellArea ! area sum - workSum(3,iRegion) = workSum(3,iRegion) + cellVolume ! volume sum - do iDataField=4,nDefinedDataFields - workSum(iDataField,iRegion) = workSum(iDataField,iRegion) + cellVolume*workArray(iDataField,iCell) ! volume-weighted sum - workMin(iDataField,iRegion) = min(workSum(iDataField,iRegion), workMin(iDataField,iRegion)) - workMax(iDataField,iRegion) = max(workSum(iDataField,iRegion), workMax(iDataField,iRegion)) - enddo + ! pack the buffer for maxval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = maxLayer(iField,iLevel,iRegion) + enddo + enddo enddo - enddo + !$omp end do + call mpas_dmpar_max_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + maxLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do - end subroutine compute_statistics + ! deallocate buffer + !$omp single + deallocate(msgBuffer) + deallocate(msgBufferReduced) + !$omp end single + + ! find min/max with region volume + !$omp do collapse(2) private(iRegion, iField) + do iRegion=1,nRegions + do iField=1,nFields + minField(iField,iRegion) = minval(minLayer(iField,:,iRegion)) + maxField(iField,iRegion) = maxval(maxLayer(iField,:,iRegion)) + enddo + enddo + !$omp end do end subroutine ocn_compute_layer_volume_weighted_averages!}}} @@ -640,6 +730,226 @@ subroutine ocn_finalize_layer_volume_weighted_averages(domain, err)!{{{ end subroutine ocn_finalize_layer_volume_weighted_averages!}}} +!*********************************************************************** +! +! routine ocn_LayerVolWgtAvgRegionMask +! +!> \brief Computes mask to define local regions +!> \author Todd Ringler, rev. Phil Jones +!> \date April 19, 2017 +!> \details +!> This is a local routine for this analysis member that computes +!> a multiplicative mask that is 1 if the point is a member of the +!> desired region and zero otherwise. +!> NOTE: This routine should be temporary and should be replaced +!> once an MPAS-wide ocean region and mask capability +!> are introduced. +! +!----------------------------------------------------------------------- + + subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & + lonCell, latCell)!{{{ + + !------------------------------------------------------------------ + ! + ! input variables + ! + !------------------------------------------------------------------ + + integer, intent(in) :: & + iBlock, &! current block in block decomposition + nCellsSolve ! number of cells owned by this block + + real(kind=RKIND), dimension(:), intent(in) :: & + lonCell, latCell ! lat, lon coordinates for each cell + + !------------------------------------------------------------------ + ! + ! output variables + ! + !------------------------------------------------------------------ + + ! output in in module variable regionMask + + !------------------------------------------------------------------ + ! + ! local variables + ! + !------------------------------------------------------------------ + + integer :: & + iCell, iRegion, &! loop indices + nCells, nRegions, nCellsMax ! array sizes + + real(kind=RKIND) :: & + dtr ! degree to radian conversion + + !------------------------------------------------------------------ + + nCellsMax = size(regionMask, dim=1) + nRegions = size(regionMask, dim=2) + + ! Currently assumes 6 regions, so fail if nRegions not large enough + if (nRegions < 6) then + !$omp master + write (stderrUnit,*) 'Abort: nRegions inconsistent ' + write (stderrUnit,*) ' : in layer volume weighted avg AM' + !$omp end master + call mpas_dmpar_global_abort(& + 'MPAS-ocean: Abort: nRegions inconsistent in layer vol wgt avg AM') + endif + + ! degrees to radians + dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND + + ! Define masks for each region + !$omp do + do iCell=1,nCellsSolve + + ! Define Arctic region (region 1) + if (latCell(iCell) < 60.0_RKIND*dtr) then + regionMask(iCell,1,iBlock) = 0.0_RKIND + else + regionMask(iCell,1,iBlock) = 1.0_RKIND + endif + + ! Equatorial region (region 2) + if (latCell(iCell) > 15.0_RKIND*dtr .or. & + latCell(iCell) < -15.0_RKIND*dtr) then + regionMask(iCell,2,iBlock) = 0.0_RKIND + else + regionMask(iCell,2,iBlock) = 1.0_RKIND + endif + + ! Southern Ocean (region 3) + if (latCell(iCell) > -50.0_RKIND*dtr) then + regionMask(iCell,3,iBlock) = 0.0_RKIND + else + regionMask(iCell,3,iBlock) = 1.0_RKIND + endif + + ! Nino 3 (region 4) + if (latCell(iCell) > 5.0_RKIND*dtr .or. & + latCell(iCell) < -5.0_RKIND*dtr .or. & + lonCell(iCell) < 210.0_RKIND*dtr .or. & + lonCell(iCell) > 270.0_RKIND*dtr) then + regionMask(iCell,4,iBlock) = 0.0_RKIND + else + regionMask(iCell,4,iBlock) = 1.0_RKIND + endif + + ! Nino 4 (region 5) + if (latCell(iCell) > 5.0_RKIND*dtr .or. & + latCell(iCell) < -5.0_RKIND*dtr .or. & + lonCell(iCell) < 160.0_RKIND*dtr .or. & + lonCell(iCell) > 210.0_RKIND*dtr) then + regionMask(iCell,5,iBlock) = 0.0_RKIND + else + regionMask(iCell,5,iBlock) = 1.0_RKIND + endif + + ! Nino 3.4 (region 6) + if (latCell(iCell) > 5.0_RKIND*dtr .or. & + latCell(iCell) < -5.0_RKIND*dtr .or. & + lonCell(iCell) < 190.0_RKIND*dtr .or. & + lonCell(iCell) > 240.0_RKIND*dtr) then + regionMask(iCell,6,iBlock) = 0.0_RKIND + else + regionMask(iCell,6,iBlock) = 1.0_RKIND + endif + + enddo + !$omp end do + + ! since number of cells varies per block, take care of any + ! extra cells in this block + + !$omp do collapse(2) + do iRegion=1,nRegions + do iCell=nCellsSolve+1,nCellsMax + regionMask(iCell,iRegion,iBlock) = 0.0_RKIND + end do + end do + !$omp end do + + !-------------------------------------------------------------------- + + end subroutine ocn_LayerVolWgtAvgRegionMask !}}} + +!*********************************************************************** +! +! routine ocn_LayerVolWgtAvgVertMask +! +!> \brief Computes mask to define active layers +!> \author Todd Ringler, rev. Phil Jones +!> \date April 19, 2017 +!> \details +!> This is a local routine for this analysis member that computes +!> a multiplicative mask that is 1 if a vertical layer is active +!> and is zero otherwise. +!> NOTE: This routine should be temporary. An MPAS ocean-wide +!> vertical mask will be used eventually. +! +!----------------------------------------------------------------------- + + subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & + vertMask)!{{{ + + !------------------------------------------------------------------ + ! + ! input variables + ! + !------------------------------------------------------------------ + + integer, intent(in) :: & + nCellsSolve ! number of cells owned by this subdomain + + integer, intent(in), dimension(:) :: & + maxLevelCell ! last active cell in each column + + !------------------------------------------------------------------ + ! + ! output variables + ! + !------------------------------------------------------------------ + + real(kind=RKIND), dimension(:,:), intent(out) :: & + vertMask ! multiplicative mask for active points + + !------------------------------------------------------------------ + ! + ! local variables + ! + !------------------------------------------------------------------ + + integer :: & + iLevel, iCell, &! loop indices + maxLevel, &! last active level in each column + nVertLevels, &! total number of vertical levels + nCells ! total number of cells + + !------------------------------------------------------------------ + + nVertLevels = size(vertMask,dim=2) + + ! compute mask assuming all levels are active down to maxLevel + + !$omp do private(iLevel,iCell,maxLevel) + do iCell=1,nCellsSolve + maxLevel = maxLevelCell(iCell) + do iLevel=1,maxLevel + vertMask(iCell,iLevel) = 1.0_RKIND + end do + do iLevel=maxLevel+1,nVertLevels + vertMask(iCell,iLevel) = 0.0_RKIND + end do + end do + !$omp end do + + !-------------------------------------------------------------------- + + end subroutine ocn_LayerVolWgtAvgVertMask !}}} + end module ocn_layer_volume_weighted_averages ! vim: foldmethod=marker From a7d08af6246c16a75bbcb52c513a7024ce77b025 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Tue, 2 May 2017 15:49:58 -0600 Subject: [PATCH 07/11] Phil performance version 2: in place --- .../mpas_ocn_layer_volume_weighted_averages.F | 352 +++++++++++++----- 1 file changed, 250 insertions(+), 102 deletions(-) diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index a977f97eb5..b932e63a38 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -253,7 +253,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pointers to data in mesh pool integer, pointer :: nVertLevels, nCells, nCellsSolve, & - nFields, nRegions + nFields, nOceanRegionsTmp integer, pointer :: index_temperature, index_salinity integer, dimension(:), pointer :: maxLevelCell real (kind=RKIND), dimension(:), pointer :: areaCell @@ -261,7 +261,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! local variables integer :: iCell, iLevel, iRegion, iTracer, iField, iBlock - real (kind=RKIND) :: maskTmp, normFactor ! for normalizing averages + real (kind=RKIND) :: normFactor ! for normalizing averages real (kind=RKIND), dimension(:,:), allocatable :: & vertMask, &!mask for active vertical levels @@ -284,7 +284,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! find the number of regions, fields and vertical levels call mpas_pool_get_dimension(domain % blocklist % dimensions, & - 'nOceanRegionsTmp', nRegions) + 'nOceanRegionsTmp', nOceanRegionsTmp) call mpas_pool_get_dimension(domain % blocklist % dimensions, & 'nLayerVolWeightedAvgFields', nFields) call mpas_pool_get_dimension(domain % blocklist % dimensions, & @@ -342,6 +342,10 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ nCellsSolve) call mpas_pool_get_dimension(block % dimensions, 'nCells', & nCells) + call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', & + nVertLevels) + call mpas_pool_get_dimension(block % dimensions, 'nOceanRegionsTmp', & + nOceanRegionsTmp) call mpas_pool_get_dimension(tracersPool, 'index_temperature', & index_temperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', & @@ -376,99 +380,141 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate arrays to hold local block masks, and fields - allocate(fieldTmp(nCellsSolve,nFields)) !$omp single - allocate(vertMask(nCells,nVertLevels)) + allocate(fieldTmp(nVertLevels,nCellsSolve)) + allocate(vertMask(nVertLevels,nCells)) !$omp end single call ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & vertMask) - !$omp do private(iLevel,iCell,iField,iRegion, & - !$omp maskTmp, fieldTmp) - do iLevel=1,nVertLevels + iField = 1 ! mask + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, vertMask, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - do iCell=1,nCellsSolve - fieldTmp(iCell,1) = vertMask (iCell,iLevel) - fieldTmp(iCell,2) = areaCell (iCell) - fieldTmp(iCell,3) = layerThickness (iLevel,iCell) - fieldTmp(iCell,4) = density (iLevel,iCell) - fieldTmp(iCell,5) = potentialDensity (iLevel,iCell) - fieldTmp(iCell,6) = BruntVaisalaFreqTop (iLevel,iCell) - fieldTmp(iCell,7) = velocityZonal (iLevel,iCell) - fieldTmp(iCell,8) = velocityMeridional (iLevel,iCell) - fieldTmp(iCell,9) = vertVelocityTop (iLevel,iCell) - fieldTmp(iCell,10) = & - activeTracers(index_temperature,iLevel,iCell) - fieldTmp(iCell,11) = & - activeTracers(index_salinity,iLevel,iCell) - fieldTmp(iCell,12) = kineticEnergyCell (iLevel,iCell) - fieldTmp(iCell,13) = relativeVorticityCell(iLevel,iCell) - fieldTmp(iCell,14) = divergence (iLevel,iCell) - fieldTmp(iCell,15) = relativeVorticityCell(iLevel,iCell)* & - relativeVorticityCell(iLevel,iCell) - end do - - ! compute min and max - do iRegion=1,nRegions - do iField=1,nFields - do iCell=1,nCellsSolve - maskTmp = vertMask(iCell,iLevel)* & - regionMask(iCell,iRegion,iBlock) - if (maskTmp /= 0.0_RKIND) then - minLayer(iField,iLevel,iRegion) = min( & - minLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) - maxLayer(iField,iLevel,iRegion) = max( & - maxLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) - endif - end do - end do - end do - - ! for the latter fields, we want volume-weighted sums - ! so replace thickness with cell volume - do iCell=1,nCellsSolve - fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) - end do - - ! compute layer averages for each region - do iRegion=1,nRegions - - ! replace field 1 with combined vertical, region mask - do iCell=1,nCellsSolve - fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & - regionMask(iCell,iRegion,iBlock) - avgLayer(1,iLevel,iRegion) = & - avgLayer(1,iLevel,iRegion) + fieldTmp(iCell,1) - end do - - ! cell area - do iCell=1,nCellsSolve - avgLayer(2,iLevel,iRegion) = & - avgLayer(2,iLevel,iRegion) + fieldTmp(iCell,2)* & - fieldTmp(iCell,1) - end do - - ! for field three, we want cell volume instead of thickness - do iCell=1,nCellsSolve - avgLayer(3,iLevel,iRegion) = & - avgLayer(3,iLevel,iRegion) + fieldTmp(iCell,3)*& - fieldTmp(iCell,1) - end do - - ! for all remaining fields, we want a volume weighted avg - do iField=4,nFields - do iCell=1,nCellsSolve - avgLayer(iField,iLevel,iRegion) = & - avgLayer(iField,iLevel,iRegion) + & - fieldTmp(iCell,iField)* & - fieldTmp(iCell,3) * & - fieldTmp(iCell,1) - end do - end do - end do ! iRegion - - end do !iLevel + iField = 2 ! cell area + !$omp do collapse(2) private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + fieldTmp(iLevel,iCell) = areaCell(iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + ! for field 3, min, max refers to layer thickness, but the + ! average is for cell volume, so include area in mask + ! (will be doing something similar for remaining fields) + iField = 3 ! cell volume + !$omp do private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + vertMask(iLevel,iCell) = vertMask(iLevel,iCell)*areaCell(iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, layerThickness, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + ! for all remaining fields, we want volume-weighted sums + ! add layer thickness to the vertMask to get a volume mask + ! use the transposed version of layerThick in fieldTmp + !$omp do private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + vertMask(iLevel,iCell) = vertMask(iLevel,iCell)* & + layerThickness(iLevel,iCell) + end do + end do + !$omp end do + + iField = 4 ! density + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, density, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 5 ! potential density + call ocn_LayerVolWgtAvgFieldMinMaxSum( avgLayer, minLayer, & + maxLayer, potentialDensity, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 6 ! Brunt Vaisala Freq + call ocn_LayerVolWgtAvgFieldMinMaxSum( avgLayer, minLayer, & + maxLayer, BruntVaisalaFreqTop, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 7 ! zonal velocity + call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & + maxLayer, velocityZonal, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 8 ! meridional velocity + call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & + maxLayer, velocityMeridional, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 9 ! vertical velocity + call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & + maxLayer, vertVelocityTop, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 10 ! temperature + !$omp do collapse(2) private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + fieldTmp(iLevel,iCell) = & + activeTracers(index_temperature,iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 11 ! salinity + !$omp do collapse(2) private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + fieldTmp(iLevel,iCell) = & + activeTracers(index_salinity,iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 12 ! kinetic energy + call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & + maxLayer, kineticEnergyCell, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 13 ! relative vorticity + call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & + maxLayer, relativeVorticityCell, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 14 ! divergence + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, divergence, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + iField = 15 ! enstrophy + !$omp do private(iCell,iLevel) + do iCell=1,nCellsSolve + do iLevel=1,nVertLevels + fieldTmp(iLevel,iCell) = relativeVorticityCell(iLevel,iCell)* & + relativeVorticityCell(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & + nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) ! deallocate fields from this block !$omp single @@ -480,7 +526,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate buffer for message passing bufferAddr=0 - bufferLength=nRegions*nFields*nVertLevels + bufferLength=nOceanRegionsTmp*nFields*nVertLevels !$omp single allocate(msgBuffer (bufferLength)) allocate(msgBufferReduced(bufferLength)) @@ -489,7 +535,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! communication for sums ! pack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -506,7 +552,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -519,7 +565,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! compute vertical sum before layer-by-layer normalization !$omp do private(iLevel,iField,iRegion,normFactor) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iField=1,nFields avgField(iField,iRegion) = 0.0_RKIND do iLevel=1,nVertLevels @@ -542,7 +588,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! normalize averages layer-by-layer !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels ! normalize all field by total volume in each layer normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) @@ -565,7 +611,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pack the buffer for minval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -580,7 +626,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -593,7 +639,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pack the buffer for maxval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -607,7 +653,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & @@ -626,7 +672,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! find min/max with region volume !$omp do collapse(2) private(iRegion, iField) - do iRegion=1,nRegions + do iRegion=1,nOceanRegionsTmp do iField=1,nFields minField(iField,iRegion) = minval(minLayer(iField,:,iRegion)) maxField(iField,iRegion) = maxval(maxLayer(iField,:,iRegion)) @@ -730,6 +776,108 @@ subroutine ocn_finalize_layer_volume_weighted_averages(domain, err)!{{{ end subroutine ocn_finalize_layer_volume_weighted_averages!}}} +!*********************************************************************** +! +! routine ocn_LayerVolWgtAvgFieldMinMaxSum +! +!> \brief Computes min, max and sum of an input field +!> \author Phil Jones +!> \date April 19, 2017 +!> \details +!> This is a local routine for this analysis member that computes +!> the min, max, and sum of a given field under the influence of a +!> mask that defines active layers and regions. The results are packed +!> into a buffer for a later reduction call. The calling routine +!> supplies the field, a precomputed region mask and a base address +!> for the buffer. +! +!----------------------------------------------------------------------- + + subroutine ocn_LayerVolWgtAvgFieldMinMaxSum( & + localSum, localMin, localMax, field, vertMask, regMask, & + nCellsSolve, numLevels, numRegions, iField, iBlock)!{{{ + + !-------------------------------------------------------------------- + ! + ! input variables + ! + !-------------------------------------------------------------------- + + integer, intent(in) :: & + numLevels, &! number of vertical levels + nCellsSolve, &! number of cells owned by this field + numRegions, &! number of regions for diagnostics + iField, &! index of field in local arrays + iBlock ! index of current block for region mask + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertMask, &! multiplicative mask defining active levels + field ! field for which statistics desired in + ! transposed order (nCells, nVertLevels) + + real (kind=RKIND), dimension(:,:,:), intent(in) :: & + regMask ! mask defining active regions + + !-------------------------------------------------------------------- + ! + ! input/output variables + ! + !-------------------------------------------------------------------- + + ! results are accumulated to account for multiple blocks + ! and are assumed to be initialized elsewhere + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + localSum, localMin, localMax ! accumulated local results + + !-------------------------------------------------------------------- + ! + ! local variables + ! + !-------------------------------------------------------------------- + + integer :: iCell, k, iRegion, iBuff ! loop indices and addresses + + real (kind=RKIND) :: maskTmp ! temp for holding mask info + + ! end preamble + !-------------------------------------------------------------------- + ! begin routine + + ! compute min, max, sum over a masked horizontal region + + !$omp do collapse(2) private(iCell,k,iRegion,maskTmp) + do iRegion=1,numRegions + do k=1,numLevels + + ! compute sum + do iCell=1,nCellsSolve + localSum(iField,k,iRegion) = & + localSum(iField,k,iRegion) + field(k,iCell)* & + vertMask(k,iCell)* & + regMask(iCell,iRegion,iBlock) + end do + + ! compute min and max + do iCell=1,nCellsSolve + maskTmp = vertMask(k,iCell)* & + regMask(iCell,iRegion,iBlock) + if (maskTmp /= 0.0_RKIND) then + localMin(iField,k,iRegion) = min( & + localMin(iField,k,iRegion),field(k,iCell)) + localMax(iField,k,iRegion) = max( & + localMax(iField,k,iRegion),field(k,iCell)) + endif + end do + + end do + end do + !$omp end do + + !-------------------------------------------------------------------- + + end subroutine ocn_LayerVolWgtAvgFieldMinMaxSum !}}} + !*********************************************************************** ! ! routine ocn_LayerVolWgtAvgRegionMask @@ -930,7 +1078,7 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & !------------------------------------------------------------------ - nVertLevels = size(vertMask,dim=2) + nVertLevels = size(vertMask,dim=1) ! compute mask assuming all levels are active down to maxLevel @@ -938,10 +1086,10 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & do iCell=1,nCellsSolve maxLevel = maxLevelCell(iCell) do iLevel=1,maxLevel - vertMask(iCell,iLevel) = 1.0_RKIND + vertMask(iLevel,iCell) = 1.0_RKIND end do do iLevel=maxLevel+1,nVertLevels - vertMask(iCell,iLevel) = 0.0_RKIND + vertMask(iLevel,iCell) = 0.0_RKIND end do end do !$omp end do From c7aabaffd0db4760bf820712702cc940d9f6e628 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Tue, 2 May 2017 15:51:47 -0600 Subject: [PATCH 08/11] Phil performance version 3: transpose --- .../mpas_ocn_layer_volume_weighted_averages.F | 335 ++++++++++++------ 1 file changed, 231 insertions(+), 104 deletions(-) diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index b932e63a38..144f207141 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -232,11 +232,15 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: tracersPool - ! pointers to result arrays + ! arrays to return results real (kind=RKIND), dimension(:,:,:), pointer :: & - avgLayer, minLayer, maxLayer + minValueWithinOceanLayerRegion, & + maxValueWithinOceanLayerRegion, & + avgValueWithinOceanLayerRegion real (kind=RKIND), dimension(:,:), pointer :: & - avgField, minField, maxField + minValueWithinOceanVolumeRegion, & + maxValueWithinOceanVolumeRegion, & + avgValueWithinOceanVolumeRegion ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -264,9 +268,14 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ real (kind=RKIND) :: normFactor ! for normalizing averages real (kind=RKIND), dimension(:,:), allocatable :: & + avgField, minField, maxField, &!results in local ordering vertMask, &!mask for active vertical levels fieldTmp !field array in local order + real (kind=RKIND), dimension(:,:,:), allocatable :: & + avgLayer, minLayer, maxLayer !results in local ordering + + ! buffers data for message passaging integer :: bufferAddr, bufferLength real (kind=RKIND), dimension(:), allocatable :: & @@ -296,27 +305,33 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ layerVolumeWeightedAverageAMPool) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'minValueWithinOceanLayerRegion', & - minLayer) + minValueWithinOceanLayerRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'maxValueWithinOceanLayerRegion', & - maxLayer) + maxValueWithinOceanLayerRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'avgValueWithinOceanLayerRegion', & - avgLayer) + avgValueWithinOceanLayerRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'minValueWithinOceanVolumeRegion', & - minField) + minValueWithinOceanVolumeRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'maxValueWithinOceanVolumeRegion', & - maxField) + maxValueWithinOceanVolumeRegion) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'avgValueWithinOceanVolumeRegion', & - avgField) - - ! initialize so this routine works for multiple blocks + avgValueWithinOceanVolumeRegion) + + ! allocate space to hold results and initialize so + ! this routine works for multiple blocks + !$omp do single + allocate(avgLayer(nVertLevels, nOceanRegionsTmp, nFields)) + allocate(minLayer(nVertLevels, nOceanRegionsTmp, nFields)) + allocate(maxLayer(nVertLevels, nOceanRegionsTmp, nFields)) avgLayer = 0.0_RKIND minLayer = 1.e20_RKIND maxLayer = -1.e20_RKIND + !$omp end do single ! loop over blocks iBlock = 0 @@ -380,9 +395,10 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate arrays to hold local block masks, and fields + ! in transposed order !$omp single - allocate(fieldTmp(nVertLevels,nCellsSolve)) - allocate(vertMask(nVertLevels,nCells)) + allocate(fieldTmp (nCellsSolve,nVertLevels)) + allocate(vertMask (nCells,nVertLevels)) !$omp end single call ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & @@ -395,9 +411,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ iField = 2 ! cell area !$omp do collapse(2) private(iCell,iLevel) - do iCell=1,nCellsSolve do iLevel=1,nVertLevels - fieldTmp(iLevel,iCell) = areaCell(iCell) + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = areaCell(iCell) end do end do !$omp end do @@ -410,63 +426,106 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! (will be doing something similar for remaining fields) iField = 3 ! cell volume !$omp do private(iCell,iLevel) - do iCell=1,nCellsSolve do iLevel=1,nVertLevels - vertMask(iLevel,iCell) = vertMask(iLevel,iCell)*areaCell(iCell) + do iCell=1,nCellsSolve + vertMask(iCell,iLevel) = vertMask(iCell,iLevel)*areaCell(iCell) + fieldTmp(iCell,iLevel) = layerThickness(iLevel,iCell) end do end do !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, layerThickness, vertMask, regionMask, & + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) ! for all remaining fields, we want volume-weighted sums ! add layer thickness to the vertMask to get a volume mask ! use the transposed version of layerThick in fieldTmp !$omp do private(iCell,iLevel) - do iCell=1,nCellsSolve do iLevel=1,nVertLevels - vertMask(iLevel,iCell) = vertMask(iLevel,iCell)* & - layerThickness(iLevel,iCell) + do iCell=1,nCellsSolve + vertMask(iCell,iLevel) = vertMask(iCell,iLevel)* & + fieldTmp(iCell,iLevel) end do end do !$omp end do iField = 4 ! density - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, density, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = density(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 5 ! potential density - call ocn_LayerVolWgtAvgFieldMinMaxSum( avgLayer, minLayer, & - maxLayer, potentialDensity, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = potentialDensity(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 6 ! Brunt Vaisala Freq - call ocn_LayerVolWgtAvgFieldMinMaxSum( avgLayer, minLayer, & - maxLayer, BruntVaisalaFreqTop, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = BruntVaisalaFreqTop(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 7 ! zonal velocity - call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & - maxLayer, velocityZonal, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = velocityZonal(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 8 ! meridional velocity - call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & - maxLayer, velocityMeridional, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = velocityMeridional(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 9 ! vertical velocity - call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & - maxLayer, vertVelocityTop, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = vertVelocityTop(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 10 ! temperature !$omp do collapse(2) private(iCell,iLevel) - do iCell=1,nCellsSolve do iLevel=1,nVertLevels - fieldTmp(iLevel,iCell) = & + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = & activeTracers(index_temperature,iLevel,iCell) end do end do @@ -477,9 +536,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ iField = 11 ! salinity !$omp do collapse(2) private(iCell,iLevel) - do iCell=1,nCellsSolve do iLevel=1,nVertLevels - fieldTmp(iLevel,iCell) = & + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = & activeTracers(index_salinity,iLevel,iCell) end do end do @@ -489,25 +548,46 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 12 ! kinetic energy - call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & - maxLayer, kineticEnergyCell, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = kineticEnergyCell(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 13 ! relative vorticity - call ocn_LayerVolWgtAvgFieldMinMaxSum(avgLayer, minLayer, & - maxLayer, relativeVorticityCell, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = relativeVorticityCell(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 14 ! divergence - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, divergence, vertMask, regionMask, & + !$omp do collapse(2) private(iCell,iLevel) + do iLevel=1,nVertLevels + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = divergence(iLevel,iCell) + end do + end do + !$omp end do + call ocn_LayerVolWgtAvgFieldMinMaxSum( & + avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) iField = 15 ! enstrophy - !$omp do private(iCell,iLevel) - do iCell=1,nCellsSolve + !$omp do collapse(2) private(iCell,iLevel) do iLevel=1,nVertLevels - fieldTmp(iLevel,iCell) = relativeVorticityCell(iLevel,iCell)* & + do iCell=1,nCellsSolve + fieldTmp(iCell,iLevel) = relativeVorticityCell(iLevel,iCell)* & relativeVorticityCell(iLevel,iCell) end do end do @@ -527,20 +607,18 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate buffer for message passing bufferAddr=0 bufferLength=nOceanRegionsTmp*nFields*nVertLevels - !$omp single allocate(msgBuffer (bufferLength)) allocate(msgBufferReduced(bufferLength)) - !$omp end single ! communication for sums ! pack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = avgLayer(iField,iLevel,iRegion) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + msgBuffer(bufferAddr) = avgLayer(iLevel,iRegion,iField) enddo enddo enddo @@ -552,37 +630,42 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - avgLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + avgLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) enddo enddo enddo !$omp end do + ! allocate space to hold results vertical means + !$omp single + allocate(avgField(nOceanRegionsTmp,nFields)) + !$omp end single + ! compute vertical sum before layer-by-layer normalization !$omp do private(iLevel,iField,iRegion,normFactor) do iRegion=1,nOceanRegionsTmp do iField=1,nFields - avgField(iField,iRegion) = 0.0_RKIND + avgField(iRegion,iField) = 0.0_RKIND do iLevel=1,nVertLevels - avgField(iField,iRegion) = avgField(iField,iRegion) & - + avgLayer(iField,iLevel,iRegion) + avgField(iRegion,iField) = avgField(iRegion,iField) & + + avgLayer(iLevel,iRegion,iField) enddo enddo - normFactor = 1.0_RKIND / max(avgField(3,iRegion),1.0e-8_RKIND) + normFactor = 1.0_RKIND / max(avgField(iRegion,3),1.0e-8_RKIND) do iField=4,nFields - avgField(iField,iRegion) = avgField(iField,iRegion)*normFactor + avgField(iRegion,iField) = avgField(iRegion,iField)*normFactor enddo ! normalize total region volume by total volume cell area - avgField(3,iRegion) = avgField(3,iRegion) & - / max(avgField(2,iRegion),1.0e-8_RKIND) + avgField(iRegion,3) = avgField(iRegion,3) & + / max(avgField(iRegion,2),1.0e-8_RKIND) ! normalize total volume cell area by total number of cells - avgField(2,iRegion) = avgField(2,iRegion) & - / max(avgField(1,iRegion),1.0e-8_RKIND) + avgField(iRegion,2) = avgField(iRegion,2) & + / max(avgField(iRegion,1),1.0e-8_RKIND) enddo !$omp end do @@ -591,32 +674,46 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels ! normalize all field by total volume in each layer - normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) + normFactor = 1.0_RKIND/max(avgLayer(iLevel,iRegion,3),1.0e-8_RKIND) do iField=4,nFields - avgLayer(iField,iLevel,iRegion) = & - avgLayer(iField,iLevel,iRegion) * normFactor + avgLayer(iLevel,iRegion,iField) = & + avgLayer(iLevel,iRegion,iField) * normFactor enddo ! normalize total layer volume by layer area - avgLayer(3,iLevel,iRegion) = avgLayer(3,iLevel,iRegion) & - / max(avgLayer(2,iLevel,iRegion),1.0e-8_RKIND) + avgLayer(iLevel,iRegion,3) = avgLayer(iLevel,iRegion,3) & + / max(avgLayer(iLevel,iRegion,2),1.0e-8_RKIND) ! normalize total layer area by number of cells in region - avgLayer(2,iLevel,iRegion) = avgLayer(2,iLevel,iRegion) & - / max(avgLayer(1,iLevel,iRegion),1.0e-8_RKIND) + avgLayer(iLevel,iRegion,2) = avgLayer(iLevel,iRegion,2) & + / max(avgLayer(iLevel,iRegion,1),1.0e-8_RKIND) enddo enddo !$omp end do + ! copy results into final location + !$omp do collapse(2) private(iLevel, iRegion, iField) + do iField=1,nFields + do iRegion=1,nOceanRegionsTmp + avgValueWithinOceanVolumeRegion(iField,iRegion) = & + avgField(iRegion,iField) + do iLevel=1,nVertLevels + avgValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & + avgLayer(iLevel,iRegion,iField) + end do + end do + end do + !$omp end do + ! now compute min/max for each field, region ! find global min/max ! pack the buffer for minval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = minLayer(iField,iLevel,iRegion) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + msgBuffer(bufferAddr) = minLayer(iLevel,iRegion,iField) enddo enddo enddo @@ -626,12 +723,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - minLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + minLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) enddo enddo enddo @@ -639,12 +736,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pack the buffer for maxval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = maxLayer(iField,iLevel,iRegion) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + msgBuffer(bufferAddr) = maxLayer(iLevel,iRegion,iField) enddo enddo enddo @@ -653,33 +750,63 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iField=1,nFields do iRegion=1,nOceanRegionsTmp do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - maxLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & + (iRegion -1)*nVertLevels + iLevel + maxLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) enddo enddo enddo !$omp end do - ! deallocate buffer + ! allocate space to hold results and free up avg fields !$omp single deallocate(msgBuffer) deallocate(msgBufferReduced) + deallocate(avgLayer) + deallocate(avgField) + allocate(minField(nOceanRegionsTmp,nFields)) + allocate(maxField(nOceanRegionsTmp,nFields)) !$omp end single ! find min/max with region volume !$omp do collapse(2) private(iRegion, iField) - do iRegion=1,nOceanRegionsTmp do iField=1,nFields - minField(iField,iRegion) = minval(minLayer(iField,:,iRegion)) - maxField(iField,iRegion) = maxval(maxLayer(iField,:,iRegion)) + do iRegion=1,nOceanRegionsTmp + minField(iRegion,iField) = minval(minLayer(:,iRegion,iField)) + maxField(iRegion,iField) = maxval(maxLayer(:,iRegion,iField)) enddo enddo !$omp end do + ! now store in final resting place + !$omp do collapse(2) private(iLevel, iRegion, iField) + do iField=1,nFields + do iRegion=1,nOceanRegionsTmp + minValueWithinOceanVolumeRegion(iField,iRegion) = & + minField(iRegion,iField) + maxValueWithinOceanVolumeRegion(iField,iRegion) = & + maxField(iRegion,iField) + do iLevel=1,nVertLevels + minValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & + minLayer(iLevel,iRegion,iField) + maxValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & + maxLayer(iLevel,iRegion,iField) + end do + end do + end do + !$omp end do + + ! deallocate remaining arrays + !$omp single + deallocate(minLayer) + deallocate(maxLayer) + deallocate(minField) + deallocate(maxField) + !$omp end single + end subroutine ocn_compute_layer_volume_weighted_averages!}}} !*********************************************************************** @@ -852,21 +979,21 @@ subroutine ocn_LayerVolWgtAvgFieldMinMaxSum( & ! compute sum do iCell=1,nCellsSolve - localSum(iField,k,iRegion) = & - localSum(iField,k,iRegion) + field(k,iCell)* & - vertMask(k,iCell)* & + localSum(k,iRegion,iField) = & + localSum(k,iRegion,iField) + field(iCell,k)* & + vertMask(iCell,k)* & regMask(iCell,iRegion,iBlock) end do ! compute min and max do iCell=1,nCellsSolve - maskTmp = vertMask(k,iCell)* & + maskTmp = vertMask(iCell,k)* & regMask(iCell,iRegion,iBlock) if (maskTmp /= 0.0_RKIND) then - localMin(iField,k,iRegion) = min( & - localMin(iField,k,iRegion),field(k,iCell)) - localMax(iField,k,iRegion) = max( & - localMax(iField,k,iRegion),field(k,iCell)) + localMin(k,iRegion,iField) = min( & + localMin(k,iRegion,iField),field(iCell,k)) + localMax(k,iRegion,iField) = max( & + localMax(k,iRegion,iField),field(iCell,k)) endif end do @@ -1078,7 +1205,7 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & !------------------------------------------------------------------ - nVertLevels = size(vertMask,dim=1) + nVertLevels = size(vertMask,dim=2) ! compute mask assuming all levels are active down to maxLevel @@ -1086,10 +1213,10 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & do iCell=1,nCellsSolve maxLevel = maxLevelCell(iCell) do iLevel=1,maxLevel - vertMask(iLevel,iCell) = 1.0_RKIND + vertMask(iCell,iLevel) = 1.0_RKIND end do do iLevel=maxLevel+1,nVertLevels - vertMask(iLevel,iCell) = 0.0_RKIND + vertMask(iCell,iLevel) = 0.0_RKIND end do end do !$omp end do From d97605e00a392e6522ac631819a2555eb4bfa3ca Mon Sep 17 00:00:00 2001 From: Anne S Berres Date: Thu, 11 May 2017 14:34:25 -0600 Subject: [PATCH 09/11] Back to Phil performance version 1: column ordering --- .../mpas_ocn_layer_volume_weighted_averages.F | 585 +++++------------- 1 file changed, 155 insertions(+), 430 deletions(-) diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index 144f207141..a977f97eb5 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -232,15 +232,11 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: tracersPool - ! arrays to return results + ! pointers to result arrays real (kind=RKIND), dimension(:,:,:), pointer :: & - minValueWithinOceanLayerRegion, & - maxValueWithinOceanLayerRegion, & - avgValueWithinOceanLayerRegion + avgLayer, minLayer, maxLayer real (kind=RKIND), dimension(:,:), pointer :: & - minValueWithinOceanVolumeRegion, & - maxValueWithinOceanVolumeRegion, & - avgValueWithinOceanVolumeRegion + avgField, minField, maxField ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -257,7 +253,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pointers to data in mesh pool integer, pointer :: nVertLevels, nCells, nCellsSolve, & - nFields, nOceanRegionsTmp + nFields, nRegions integer, pointer :: index_temperature, index_salinity integer, dimension(:), pointer :: maxLevelCell real (kind=RKIND), dimension(:), pointer :: areaCell @@ -265,17 +261,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! local variables integer :: iCell, iLevel, iRegion, iTracer, iField, iBlock - real (kind=RKIND) :: normFactor ! for normalizing averages + real (kind=RKIND) :: maskTmp, normFactor ! for normalizing averages real (kind=RKIND), dimension(:,:), allocatable :: & - avgField, minField, maxField, &!results in local ordering vertMask, &!mask for active vertical levels fieldTmp !field array in local order - real (kind=RKIND), dimension(:,:,:), allocatable :: & - avgLayer, minLayer, maxLayer !results in local ordering - - ! buffers data for message passaging integer :: bufferAddr, bufferLength real (kind=RKIND), dimension(:), allocatable :: & @@ -293,7 +284,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! find the number of regions, fields and vertical levels call mpas_pool_get_dimension(domain % blocklist % dimensions, & - 'nOceanRegionsTmp', nOceanRegionsTmp) + 'nOceanRegionsTmp', nRegions) call mpas_pool_get_dimension(domain % blocklist % dimensions, & 'nLayerVolWeightedAvgFields', nFields) call mpas_pool_get_dimension(domain % blocklist % dimensions, & @@ -305,33 +296,27 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ layerVolumeWeightedAverageAMPool) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'minValueWithinOceanLayerRegion', & - minValueWithinOceanLayerRegion) + minLayer) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'maxValueWithinOceanLayerRegion', & - maxValueWithinOceanLayerRegion) + maxLayer) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'avgValueWithinOceanLayerRegion', & - avgValueWithinOceanLayerRegion) + avgLayer) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'minValueWithinOceanVolumeRegion', & - minValueWithinOceanVolumeRegion) + minField) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'maxValueWithinOceanVolumeRegion', & - maxValueWithinOceanVolumeRegion) + maxField) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & 'avgValueWithinOceanVolumeRegion', & - avgValueWithinOceanVolumeRegion) - - ! allocate space to hold results and initialize so - ! this routine works for multiple blocks - !$omp do single - allocate(avgLayer(nVertLevels, nOceanRegionsTmp, nFields)) - allocate(minLayer(nVertLevels, nOceanRegionsTmp, nFields)) - allocate(maxLayer(nVertLevels, nOceanRegionsTmp, nFields)) + avgField) + + ! initialize so this routine works for multiple blocks avgLayer = 0.0_RKIND minLayer = 1.e20_RKIND maxLayer = -1.e20_RKIND - !$omp end do single ! loop over blocks iBlock = 0 @@ -357,10 +342,6 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ nCellsSolve) call mpas_pool_get_dimension(block % dimensions, 'nCells', & nCells) - call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', & - nVertLevels) - call mpas_pool_get_dimension(block % dimensions, 'nOceanRegionsTmp', & - nOceanRegionsTmp) call mpas_pool_get_dimension(tracersPool, 'index_temperature', & index_temperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', & @@ -395,206 +376,99 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate arrays to hold local block masks, and fields - ! in transposed order + allocate(fieldTmp(nCellsSolve,nFields)) !$omp single - allocate(fieldTmp (nCellsSolve,nVertLevels)) - allocate(vertMask (nCells,nVertLevels)) + allocate(vertMask(nCells,nVertLevels)) !$omp end single call ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & vertMask) - iField = 1 ! mask - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, vertMask, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 2 ! cell area - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = areaCell(iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - ! for field 3, min, max refers to layer thickness, but the - ! average is for cell volume, so include area in mask - ! (will be doing something similar for remaining fields) - iField = 3 ! cell volume - !$omp do private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - vertMask(iCell,iLevel) = vertMask(iCell,iLevel)*areaCell(iCell) - fieldTmp(iCell,iLevel) = layerThickness(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - ! for all remaining fields, we want volume-weighted sums - ! add layer thickness to the vertMask to get a volume mask - ! use the transposed version of layerThick in fieldTmp - !$omp do private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - vertMask(iCell,iLevel) = vertMask(iCell,iLevel)* & - fieldTmp(iCell,iLevel) - end do - end do - !$omp end do - - iField = 4 ! density - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = density(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 5 ! potential density - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = potentialDensity(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 6 ! Brunt Vaisala Freq - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = BruntVaisalaFreqTop(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 7 ! zonal velocity - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = velocityZonal(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 8 ! meridional velocity - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = velocityMeridional(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 9 ! vertical velocity - !$omp do collapse(2) private(iCell,iLevel) + !$omp do private(iLevel,iCell,iField,iRegion, & + !$omp maskTmp, fieldTmp) do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = vertVelocityTop(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 10 ! temperature - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = & - activeTracers(index_temperature,iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 11 ! salinity - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = & - activeTracers(index_salinity,iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 12 ! kinetic energy - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = kineticEnergyCell(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 13 ! relative vorticity - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = relativeVorticityCell(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 14 ! divergence - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = divergence(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) - - iField = 15 ! enstrophy - !$omp do collapse(2) private(iCell,iLevel) - do iLevel=1,nVertLevels - do iCell=1,nCellsSolve - fieldTmp(iCell,iLevel) = relativeVorticityCell(iLevel,iCell)* & - relativeVorticityCell(iLevel,iCell) - end do - end do - !$omp end do - call ocn_LayerVolWgtAvgFieldMinMaxSum( & - avgLayer, minLayer, maxLayer, fieldTmp, vertMask, regionMask, & - nCellsSolve, nVertLevels, nOceanRegionsTmp, iField, iBlock) + + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask (iCell,iLevel) + fieldTmp(iCell,2) = areaCell (iCell) + fieldTmp(iCell,3) = layerThickness (iLevel,iCell) + fieldTmp(iCell,4) = density (iLevel,iCell) + fieldTmp(iCell,5) = potentialDensity (iLevel,iCell) + fieldTmp(iCell,6) = BruntVaisalaFreqTop (iLevel,iCell) + fieldTmp(iCell,7) = velocityZonal (iLevel,iCell) + fieldTmp(iCell,8) = velocityMeridional (iLevel,iCell) + fieldTmp(iCell,9) = vertVelocityTop (iLevel,iCell) + fieldTmp(iCell,10) = & + activeTracers(index_temperature,iLevel,iCell) + fieldTmp(iCell,11) = & + activeTracers(index_salinity,iLevel,iCell) + fieldTmp(iCell,12) = kineticEnergyCell (iLevel,iCell) + fieldTmp(iCell,13) = relativeVorticityCell(iLevel,iCell) + fieldTmp(iCell,14) = divergence (iLevel,iCell) + fieldTmp(iCell,15) = relativeVorticityCell(iLevel,iCell)* & + relativeVorticityCell(iLevel,iCell) + end do + + ! compute min and max + do iRegion=1,nRegions + do iField=1,nFields + do iCell=1,nCellsSolve + maskTmp = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + if (maskTmp /= 0.0_RKIND) then + minLayer(iField,iLevel,iRegion) = min( & + minLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + maxLayer(iField,iLevel,iRegion) = max( & + maxLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + endif + end do + end do + end do + + ! for the latter fields, we want volume-weighted sums + ! so replace thickness with cell volume + do iCell=1,nCellsSolve + fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) + end do + + ! compute layer averages for each region + do iRegion=1,nRegions + + ! replace field 1 with combined vertical, region mask + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + avgLayer(1,iLevel,iRegion) = & + avgLayer(1,iLevel,iRegion) + fieldTmp(iCell,1) + end do + + ! cell area + do iCell=1,nCellsSolve + avgLayer(2,iLevel,iRegion) = & + avgLayer(2,iLevel,iRegion) + fieldTmp(iCell,2)* & + fieldTmp(iCell,1) + end do + + ! for field three, we want cell volume instead of thickness + do iCell=1,nCellsSolve + avgLayer(3,iLevel,iRegion) = & + avgLayer(3,iLevel,iRegion) + fieldTmp(iCell,3)*& + fieldTmp(iCell,1) + end do + + ! for all remaining fields, we want a volume weighted avg + do iField=4,nFields + do iCell=1,nCellsSolve + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) + & + fieldTmp(iCell,iField)* & + fieldTmp(iCell,3) * & + fieldTmp(iCell,1) + end do + end do + end do ! iRegion + + end do !iLevel ! deallocate fields from this block !$omp single @@ -606,19 +480,21 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! allocate buffer for message passing bufferAddr=0 - bufferLength=nOceanRegionsTmp*nFields*nVertLevels + bufferLength=nRegions*nFields*nVertLevels + !$omp single allocate(msgBuffer (bufferLength)) allocate(msgBufferReduced(bufferLength)) + !$omp end single ! communication for sums ! pack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - msgBuffer(bufferAddr) = avgLayer(iLevel,iRegion,iField) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = avgLayer(iField,iLevel,iRegion) enddo enddo enddo @@ -630,90 +506,71 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - avgLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + avgLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo !$omp end do - ! allocate space to hold results vertical means - !$omp single - allocate(avgField(nOceanRegionsTmp,nFields)) - !$omp end single - ! compute vertical sum before layer-by-layer normalization !$omp do private(iLevel,iField,iRegion,normFactor) - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iField=1,nFields - avgField(iRegion,iField) = 0.0_RKIND + avgField(iField,iRegion) = 0.0_RKIND do iLevel=1,nVertLevels - avgField(iRegion,iField) = avgField(iRegion,iField) & - + avgLayer(iLevel,iRegion,iField) + avgField(iField,iRegion) = avgField(iField,iRegion) & + + avgLayer(iField,iLevel,iRegion) enddo enddo - normFactor = 1.0_RKIND / max(avgField(iRegion,3),1.0e-8_RKIND) + normFactor = 1.0_RKIND / max(avgField(3,iRegion),1.0e-8_RKIND) do iField=4,nFields - avgField(iRegion,iField) = avgField(iRegion,iField)*normFactor + avgField(iField,iRegion) = avgField(iField,iRegion)*normFactor enddo ! normalize total region volume by total volume cell area - avgField(iRegion,3) = avgField(iRegion,3) & - / max(avgField(iRegion,2),1.0e-8_RKIND) + avgField(3,iRegion) = avgField(3,iRegion) & + / max(avgField(2,iRegion),1.0e-8_RKIND) ! normalize total volume cell area by total number of cells - avgField(iRegion,2) = avgField(iRegion,2) & - / max(avgField(iRegion,1),1.0e-8_RKIND) + avgField(2,iRegion) = avgField(2,iRegion) & + / max(avgField(1,iRegion),1.0e-8_RKIND) enddo !$omp end do ! normalize averages layer-by-layer !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels ! normalize all field by total volume in each layer - normFactor = 1.0_RKIND/max(avgLayer(iLevel,iRegion,3),1.0e-8_RKIND) + normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) do iField=4,nFields - avgLayer(iLevel,iRegion,iField) = & - avgLayer(iLevel,iRegion,iField) * normFactor + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) * normFactor enddo ! normalize total layer volume by layer area - avgLayer(iLevel,iRegion,3) = avgLayer(iLevel,iRegion,3) & - / max(avgLayer(iLevel,iRegion,2),1.0e-8_RKIND) + avgLayer(3,iLevel,iRegion) = avgLayer(3,iLevel,iRegion) & + / max(avgLayer(2,iLevel,iRegion),1.0e-8_RKIND) ! normalize total layer area by number of cells in region - avgLayer(iLevel,iRegion,2) = avgLayer(iLevel,iRegion,2) & - / max(avgLayer(iLevel,iRegion,1),1.0e-8_RKIND) + avgLayer(2,iLevel,iRegion) = avgLayer(2,iLevel,iRegion) & + / max(avgLayer(1,iLevel,iRegion),1.0e-8_RKIND) enddo enddo !$omp end do - ! copy results into final location - !$omp do collapse(2) private(iLevel, iRegion, iField) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp - avgValueWithinOceanVolumeRegion(iField,iRegion) = & - avgField(iRegion,iField) - do iLevel=1,nVertLevels - avgValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & - avgLayer(iLevel,iRegion,iField) - end do - end do - end do - !$omp end do - ! now compute min/max for each field, region ! find global min/max ! pack the buffer for minval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - msgBuffer(bufferAddr) = minLayer(iLevel,iRegion,iField) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = minLayer(iField,iLevel,iRegion) enddo enddo enddo @@ -723,12 +580,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - minLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + minLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo @@ -736,12 +593,12 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pack the buffer for maxval !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - msgBuffer(bufferAddr) = maxLayer(iLevel,iRegion,iField) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = maxLayer(iField,iLevel,iRegion) enddo enddo enddo @@ -750,63 +607,33 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ msgBufferReduced ) ! unpack the buffer !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp + do iRegion=1,nRegions do iLevel=1,nVertLevels - bufferAddr = (iField -1)*nOceanRegionsTmp*nVertLevels + & - (iRegion -1)*nVertLevels + iLevel - maxLayer(iLevel,iRegion,iField)=msgBufferReduced(bufferAddr) + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + maxLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo !$omp end do - ! allocate space to hold results and free up avg fields + ! deallocate buffer !$omp single deallocate(msgBuffer) deallocate(msgBufferReduced) - deallocate(avgLayer) - deallocate(avgField) - allocate(minField(nOceanRegionsTmp,nFields)) - allocate(maxField(nOceanRegionsTmp,nFields)) !$omp end single ! find min/max with region volume !$omp do collapse(2) private(iRegion, iField) + do iRegion=1,nRegions do iField=1,nFields - do iRegion=1,nOceanRegionsTmp - minField(iRegion,iField) = minval(minLayer(:,iRegion,iField)) - maxField(iRegion,iField) = maxval(maxLayer(:,iRegion,iField)) + minField(iField,iRegion) = minval(minLayer(iField,:,iRegion)) + maxField(iField,iRegion) = maxval(maxLayer(iField,:,iRegion)) enddo enddo !$omp end do - ! now store in final resting place - !$omp do collapse(2) private(iLevel, iRegion, iField) - do iField=1,nFields - do iRegion=1,nOceanRegionsTmp - minValueWithinOceanVolumeRegion(iField,iRegion) = & - minField(iRegion,iField) - maxValueWithinOceanVolumeRegion(iField,iRegion) = & - maxField(iRegion,iField) - do iLevel=1,nVertLevels - minValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & - minLayer(iLevel,iRegion,iField) - maxValueWithinOceanLayerRegion(iField,iLevel,iRegion) = & - maxLayer(iLevel,iRegion,iField) - end do - end do - end do - !$omp end do - - ! deallocate remaining arrays - !$omp single - deallocate(minLayer) - deallocate(maxLayer) - deallocate(minField) - deallocate(maxField) - !$omp end single - end subroutine ocn_compute_layer_volume_weighted_averages!}}} !*********************************************************************** @@ -903,108 +730,6 @@ subroutine ocn_finalize_layer_volume_weighted_averages(domain, err)!{{{ end subroutine ocn_finalize_layer_volume_weighted_averages!}}} -!*********************************************************************** -! -! routine ocn_LayerVolWgtAvgFieldMinMaxSum -! -!> \brief Computes min, max and sum of an input field -!> \author Phil Jones -!> \date April 19, 2017 -!> \details -!> This is a local routine for this analysis member that computes -!> the min, max, and sum of a given field under the influence of a -!> mask that defines active layers and regions. The results are packed -!> into a buffer for a later reduction call. The calling routine -!> supplies the field, a precomputed region mask and a base address -!> for the buffer. -! -!----------------------------------------------------------------------- - - subroutine ocn_LayerVolWgtAvgFieldMinMaxSum( & - localSum, localMin, localMax, field, vertMask, regMask, & - nCellsSolve, numLevels, numRegions, iField, iBlock)!{{{ - - !-------------------------------------------------------------------- - ! - ! input variables - ! - !-------------------------------------------------------------------- - - integer, intent(in) :: & - numLevels, &! number of vertical levels - nCellsSolve, &! number of cells owned by this field - numRegions, &! number of regions for diagnostics - iField, &! index of field in local arrays - iBlock ! index of current block for region mask - - real (kind=RKIND), dimension(:,:), intent(in) :: & - vertMask, &! multiplicative mask defining active levels - field ! field for which statistics desired in - ! transposed order (nCells, nVertLevels) - - real (kind=RKIND), dimension(:,:,:), intent(in) :: & - regMask ! mask defining active regions - - !-------------------------------------------------------------------- - ! - ! input/output variables - ! - !-------------------------------------------------------------------- - - ! results are accumulated to account for multiple blocks - ! and are assumed to be initialized elsewhere - - real (kind=RKIND), dimension(:,:,:), intent(inout) :: & - localSum, localMin, localMax ! accumulated local results - - !-------------------------------------------------------------------- - ! - ! local variables - ! - !-------------------------------------------------------------------- - - integer :: iCell, k, iRegion, iBuff ! loop indices and addresses - - real (kind=RKIND) :: maskTmp ! temp for holding mask info - - ! end preamble - !-------------------------------------------------------------------- - ! begin routine - - ! compute min, max, sum over a masked horizontal region - - !$omp do collapse(2) private(iCell,k,iRegion,maskTmp) - do iRegion=1,numRegions - do k=1,numLevels - - ! compute sum - do iCell=1,nCellsSolve - localSum(k,iRegion,iField) = & - localSum(k,iRegion,iField) + field(iCell,k)* & - vertMask(iCell,k)* & - regMask(iCell,iRegion,iBlock) - end do - - ! compute min and max - do iCell=1,nCellsSolve - maskTmp = vertMask(iCell,k)* & - regMask(iCell,iRegion,iBlock) - if (maskTmp /= 0.0_RKIND) then - localMin(k,iRegion,iField) = min( & - localMin(k,iRegion,iField),field(iCell,k)) - localMax(k,iRegion,iField) = max( & - localMax(k,iRegion,iField),field(iCell,k)) - endif - end do - - end do - end do - !$omp end do - - !-------------------------------------------------------------------- - - end subroutine ocn_LayerVolWgtAvgFieldMinMaxSum !}}} - !*********************************************************************** ! ! routine ocn_LayerVolWgtAvgRegionMask From 187850de3980369fbf925a36b697b5175e03e4a3 Mon Sep 17 00:00:00 2001 From: Anne S Berres Date: Thu, 4 May 2017 14:25:41 -0600 Subject: [PATCH 10/11] Region mask support for Layer Volume Weighted Averages. Hard-wired regions and region masks can be turned on/off individually in namelist. --- ...egistry_layer_volume_weighted_averages.xml | 303 ++++++- .../mpas_ocn_layer_volume_weighted_averages.F | 830 ++++++++++++------ 2 files changed, 881 insertions(+), 252 deletions(-) diff --git a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml index 1f3aa5137b..9fe27ace2a 100644 --- a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml +++ b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml @@ -2,7 +2,7 @@ - @@ -27,6 +27,14 @@ description="Name of the string that should be tied to the layer volume weighted average analysis member" possible_values="Any existing stream name or 'none'" /> + + @@ -314,6 +322,288 @@ description="Average relative enstrophy within region volume" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index a977f97eb5..48f045c09c 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -115,7 +115,7 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ integer, pointer :: & nCells, &! number of cells in a given block nCellsSolve, &! number of cells owned by a given block - nRegions ! number of ocean regions defined + nRegionsTmp ! number of ocean regions defined real (kind=RKIND), dimension(:), pointer :: & lonCell, latCell ! latitude, longitude of each cell @@ -134,7 +134,7 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ ! retrieve number of regions call mpas_pool_get_dimension(domain % blocklist % dimensions, & - 'nOceanRegionsTmp', nRegions) + 'nOceanRegionsTmp', nRegionsTmp) ! count up the number of blocks and the max size of each domain numBlocks = 0 @@ -151,7 +151,7 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ end do ! allocate the region mask array - allocate(regionMask(maxCells,nRegions,numBlocks)) + allocate(regionMask(maxCells,nRegionsTmp,numBlocks)) ! compute the region mask for each block iBlock = 0 @@ -181,8 +181,8 @@ end subroutine ocn_init_layer_volume_weighted_averages!}}} ! routine ocn_compute_layer_volume_weighted_averages ! !> \brief Compute MPAS-Ocean analysis member -!> \author Todd Ringler -!> \date April 24, 2015 +!> \author Todd Ringler, Anne Berres +!> \date May 8, 2017 !> \details !> This routine conducts all computation required for this !> MPAS-Ocean analysis member. @@ -231,12 +231,13 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ type (mpas_pool_type), pointer :: diagnosticsPool type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: tracersPool + type (mpas_pool_type), pointer :: regionPool ! pointers to result arrays real (kind=RKIND), dimension(:,:,:), pointer :: & - avgLayer, minLayer, maxLayer + avgLayer, minLayer, maxLayer, avgLayerRF, minLayerRF, maxLayerRF real (kind=RKIND), dimension(:,:), pointer :: & - avgField, minField, maxField + avgField, minField, maxField, avgFieldRF, minFieldRF, maxFieldRF ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -253,10 +254,11 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pointers to data in mesh pool integer, pointer :: nVertLevels, nCells, nCellsSolve, & - nFields, nRegions + nFields, nRegionsTmp integer, pointer :: index_temperature, index_salinity integer, dimension(:), pointer :: maxLevelCell real (kind=RKIND), dimension(:), pointer :: areaCell + logical, pointer :: predefRegions ! local variables integer :: iCell, iLevel, iRegion, iTracer, iField, iBlock @@ -272,6 +274,15 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ real (kind=RKIND), dimension(:), allocatable :: & msgBuffer, msgBufferReduced + !!! region file variables + integer, pointer :: regionsInAddGroup + integer :: regionGroupOffset, regionGroupNumber, i, j + character (len=STRKIND), dimension(:), pointer :: regionNames, regionGroupNames + integer, dimension(:,:), pointer :: regionCellMasks, regionsInGroup + integer, dimension(:), pointer :: nRegionsInGroup + integer, pointer ::nRegions, nRegionGroups, maxRegionsInGroup + character (len=STRKIND), pointer :: additionalRegion + ! end of preamble !----------------------------------------------------------------- ! begin code @@ -282,41 +293,125 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! set highest level pointer dminfo = domain % dminfo - ! find the number of regions, fields and vertical levels - call mpas_pool_get_dimension(domain % blocklist % dimensions, & - 'nOceanRegionsTmp', nRegions) + ! find the number of fields and vertical levels call mpas_pool_get_dimension(domain % blocklist % dimensions, & 'nLayerVolWeightedAvgFields', nFields) call mpas_pool_get_dimension(domain % blocklist % dimensions, & 'nVertLevels', nVertLevels) - ! get pointers to analysis member arrays to store results + call mpas_pool_get_subpool(domain % blocklist % structs, & 'layerVolumeWeightedAverageAM', & layerVolumeWeightedAverageAMPool) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'minValueWithinOceanLayerRegion', & - minLayer) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'maxValueWithinOceanLayerRegion', & - maxLayer) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'avgValueWithinOceanLayerRegion', & - avgLayer) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'minValueWithinOceanVolumeRegion', & - minField) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'maxValueWithinOceanVolumeRegion', & - maxField) - call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'avgValueWithinOceanVolumeRegion', & - avgField) + + ! region config for LVWA + call mpas_pool_get_config(domain % configs, & + 'config_AM_layerVolumeWeightedAverage_compute_predef_regions', & + predefRegions) + call mpas_pool_get_config(domain % configs, & + 'config_AM_layerVolumeWeightedAverage_region_group', & + additionalRegion) + + if (predefRegions == .true.) then + ! get pointers to analysis member arrays to store results + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanLayerRegion', & + minLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanLayerRegion', & + maxLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanLayerRegion', & + avgLayer) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanVolumeRegion', & + minField) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanVolumeRegion', & + maxField) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanVolumeRegion', & + avgField) + + ! find the number of regions + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nOceanRegionsTmp', nRegionsTmp) + endif + + if (additionalRegion /= '') then + + ! get pointers to analysis member arrays to store results + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanLayerRF', & + minLayerRF) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanLayerRF', & + maxLayerRF) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanLayerRF', & + avgLayerRF) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanVolumeRF', & + minFieldRF) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanVolumeRF', & + maxFieldRF) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanVolumeRF', & + avgFieldRF) + + !!! region file variables + ! if a region file is given and a region is selected, set up region file version + ! region file dimensions + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nRegions', nRegions) + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nRegionGroups', nRegionGroups) + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'maxRegionsInGroup', maxRegionsInGroup) + + ! get region file dimensions + call mpas_pool_get_subpool(domain % blocklist % structs, & + 'regions', regionPool) + call mpas_pool_get_array(regionPool, 'regionsInGroup', & + regionsInGroup) + call mpas_pool_get_array(regionPool, 'nRegionsInGroup', & + nRegionsInGroup) + call mpas_pool_get_array(regionPool, 'regionNames', & + regionNames) + call mpas_pool_get_array(regionPool, 'regionGroupNames', & + regionGroupNames) + call mpas_pool_get_array(regionPool, 'regionCellMasks', & + regionCellMasks) + + regionGroupOffset = 1 + regionGroupNumber = 0 + + ! region preliminaries + ! figure out the region group number that matches the configured additional region's name + do i=1,nRegionGroups + if (regionGroupNames(i) .eq. additionalRegion) then + regionGroupNumber = i + ! determine offset to compensate for several region groups in the + ! regions file + do j=1,i-1 + regionGroupOffset = regionGroupOffset + nRegionsInGroup(j) + enddo !j + endif + enddo !i + + ! find the number of regions + nRegions => nRegionsInGroup(regionGroupNumber) + endif + ! initialize so this routine works for multiple blocks avgLayer = 0.0_RKIND minLayer = 1.e20_RKIND maxLayer = -1.e20_RKIND + avgLayerRF = 0.0_RKIND + minLayerRF = 1.e20_RKIND + maxLayerRF = -1.e20_RKIND ! loop over blocks iBlock = 0 @@ -381,6 +476,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ allocate(vertMask(nCells,nVertLevels)) !$omp end single + call ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & vertMask) @@ -393,7 +489,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ fieldTmp(iCell,2) = areaCell (iCell) fieldTmp(iCell,3) = layerThickness (iLevel,iCell) fieldTmp(iCell,4) = density (iLevel,iCell) - fieldTmp(iCell,5) = potentialDensity (iLevel,iCell) + fieldTmp(iCell,5) = potentialDensity (iLevel,iCell) fieldTmp(iCell,6) = BruntVaisalaFreqTop (iLevel,iCell) fieldTmp(iCell,7) = velocityZonal (iLevel,iCell) fieldTmp(iCell,8) = velocityMeridional (iLevel,iCell) @@ -405,69 +501,131 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ fieldTmp(iCell,12) = kineticEnergyCell (iLevel,iCell) fieldTmp(iCell,13) = relativeVorticityCell(iLevel,iCell) fieldTmp(iCell,14) = divergence (iLevel,iCell) - fieldTmp(iCell,15) = relativeVorticityCell(iLevel,iCell)* & + fieldTmp(iCell,15) = relativeVorticityCell(iLevel,iCell)* & relativeVorticityCell(iLevel,iCell) end do - - ! compute min and max - do iRegion=1,nRegions - do iField=1,nFields - do iCell=1,nCellsSolve - maskTmp = vertMask(iCell,iLevel)* & - regionMask(iCell,iRegion,iBlock) - if (maskTmp /= 0.0_RKIND) then - minLayer(iField,iLevel,iRegion) = min( & - minLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) - maxLayer(iField,iLevel,iRegion) = max( & - maxLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) - endif - end do - end do - end do - - ! for the latter fields, we want volume-weighted sums - ! so replace thickness with cell volume - do iCell=1,nCellsSolve - fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) - end do - - ! compute layer averages for each region - do iRegion=1,nRegions - ! replace field 1 with combined vertical, region mask + if(preDefRegions == .true.) then + ! compute min and max + do iRegion=1,nRegionsTmp + do iField=1,nFields do iCell=1,nCellsSolve - fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & - regionMask(iCell,iRegion,iBlock) - avgLayer(1,iLevel,iRegion) = & - avgLayer(1,iLevel,iRegion) + fieldTmp(iCell,1) + maskTmp = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + if (maskTmp /= 0.0_RKIND) then + minLayer(iField,iLevel,iRegion) = min( & + minLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + maxLayer(iField,iLevel,iRegion) = max( & + maxLayer(iField,iLevel,iRegion),fieldTmp(iCell,iField)) + endif + end do + end do end do - ! cell area + ! for the latter fields, we want volume-weighted sums + ! so replace thickness with cell volume do iCell=1,nCellsSolve - avgLayer(2,iLevel,iRegion) = & - avgLayer(2,iLevel,iRegion) + fieldTmp(iCell,2)* & - fieldTmp(iCell,1) + fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) end do - ! for field three, we want cell volume instead of thickness + ! compute layer averages for each region + do iRegion=1,nRegionsTmp + + ! replace field 1 with combined vertical, region mask + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & + regionMask(iCell,iRegion,iBlock) + avgLayer(1,iLevel,iRegion) = & + avgLayer(1,iLevel,iRegion) + fieldTmp(iCell,1) + end do + + ! cell area + do iCell=1,nCellsSolve + avgLayer(2,iLevel,iRegion) = & + avgLayer(2,iLevel,iRegion) + fieldTmp(iCell,2)* & + fieldTmp(iCell,1) + end do + + ! for field three, we want cell volume instead of thickness + do iCell=1,nCellsSolve + avgLayer(3,iLevel,iRegion) = & + avgLayer(3,iLevel,iRegion) + fieldTmp(iCell,3)*& + fieldTmp(iCell,1) + end do + + ! for all remaining fields, we want a volume weighted avg + do iField=4,nFields + do iCell=1,nCellsSolve + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) + & + fieldTmp(iCell,iField)* & + fieldTmp(iCell,3) * & + fieldTmp(iCell,1) + end do + end do + end do ! iRegion + endif ! hard-wired regions + + if(additionalRegion /= '') then + ! compute min and max + do iRegion=1,nRegions + do iField=1,nFields do iCell=1,nCellsSolve - avgLayer(3,iLevel,iRegion) = & - avgLayer(3,iLevel,iRegion) + fieldTmp(iCell,3)*& - fieldTmp(iCell,1) + maskTmp = vertMask(iCell,iLevel) * regionCellMasks(iRegion,iCell) + if (maskTmp /= 0.0_RKIND) then + minLayerRF(iField,iLevel,iRegion) = min( & + minLayerRF(iField,iLevel,iRegion), fieldTmp(iCell,iField)) + maxLayerRF(iField,iLevel,iRegion) = max( & + maxLayerRF(iField,iLevel,iRegion), fieldTmp(iCell,iField)) + endif + end do + end do end do - ! for all remaining fields, we want a volume weighted avg - do iField=4,nFields + ! for the latter fields, we want volume-weighted sums + ! so replace thickness with cell volume do iCell=1,nCellsSolve - avgLayer(iField,iLevel,iRegion) = & - avgLayer(iField,iLevel,iRegion) + & - fieldTmp(iCell,iField)* & - fieldTmp(iCell,3) * & - fieldTmp(iCell,1) - end do + fieldTmp(iCell,3) = fieldTmp(iCell,2)*fieldTmp(iCell,3) end do - end do ! iRegion - + + ! compute layer averages for each region + do iRegion=1,nRegions + + ! replace field 1 with combined vertical, region mask + do iCell=1,nCellsSolve + fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & + regionCellMasks(iRegion,iCell) + avgLayerRF(1,iLevel,iRegion) = & + avgLayerRF(1,iLevel,iRegion) + fieldTmp(iCell,1) + end do + + ! cell area + do iCell=1,nCellsSolve + avgLayerRF(2,iLevel,iRegion) = & + avgLayerRF(2,iLevel,iRegion) + fieldTmp(iCell,2)* & + fieldTmp(iCell,1) + end do + + ! for field three, we want cell volume instead of thickness + do iCell=1,nCellsSolve + avgLayerRF(3,iLevel,iRegion) = & + avgLayerRF(3,iLevel,iRegion) + fieldTmp(iCell,3)*& + fieldTmp(iCell,1) + end do + + ! for all remaining fields, we want a volume weighted avg + do iField=4,nFields + do iCell=1,nCellsSolve + avgLayerRF(iField,iLevel,iRegion) = & + avgLayerRF(iField,iLevel,iRegion) + & + fieldTmp(iCell,iField)* & + fieldTmp(iCell,3) * & + fieldTmp(iCell,1) + end do + end do + end do ! iRegion + endif ! region mask file + end do !iLevel ! deallocate fields from this block @@ -475,165 +633,332 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ deallocate(fieldTmp) deallocate(vertMask) !$omp end single + block => block % next end do - ! allocate buffer for message passing - bufferAddr=0 - bufferLength=nRegions*nFields*nVertLevels - !$omp single - allocate(msgBuffer (bufferLength)) - allocate(msgBufferReduced(bufferLength)) - !$omp end single - - ! communication for sums - ! pack the buffer - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = avgLayer(iField,iLevel,iRegion) - enddo - enddo - enddo - !$omp end do - call mpas_dmpar_sum_real_array(dminfo, bufferLength, msgBuffer, & - msgBufferReduced ) - - ! compute normalized averages - - ! unpack the buffer - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - avgLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) - enddo - enddo - enddo - !$omp end do - - ! compute vertical sum before layer-by-layer normalization - !$omp do private(iLevel,iField,iRegion,normFactor) - do iRegion=1,nRegions - do iField=1,nFields - avgField(iField,iRegion) = 0.0_RKIND - do iLevel=1,nVertLevels - avgField(iField,iRegion) = avgField(iField,iRegion) & - + avgLayer(iField,iLevel,iRegion) - enddo - enddo - normFactor = 1.0_RKIND / max(avgField(3,iRegion),1.0e-8_RKIND) - do iField=4,nFields - avgField(iField,iRegion) = avgField(iField,iRegion)*normFactor - enddo - ! normalize total region volume by total volume cell area - avgField(3,iRegion) = avgField(3,iRegion) & - / max(avgField(2,iRegion),1.0e-8_RKIND) - ! normalize total volume cell area by total number of cells - avgField(2,iRegion) = avgField(2,iRegion) & - / max(avgField(1,iRegion),1.0e-8_RKIND) - enddo - !$omp end do - - ! normalize averages layer-by-layer - !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - ! normalize all field by total volume in each layer - normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) - do iField=4,nFields - avgLayer(iField,iLevel,iRegion) = & - avgLayer(iField,iLevel,iRegion) * normFactor - enddo - ! normalize total layer volume by layer area - avgLayer(3,iLevel,iRegion) = avgLayer(3,iLevel,iRegion) & - / max(avgLayer(2,iLevel,iRegion),1.0e-8_RKIND) - ! normalize total layer area by number of cells in region - avgLayer(2,iLevel,iRegion) = avgLayer(2,iLevel,iRegion) & - / max(avgLayer(1,iLevel,iRegion),1.0e-8_RKIND) - enddo - enddo - !$omp end do - - ! now compute min/max for each field, region - ! find global min/max - - ! pack the buffer for minval - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = minLayer(iField,iLevel,iRegion) - enddo - enddo - enddo - !$omp end do - ! communicate for global min - call mpas_dmpar_min_real_array(dminfo, bufferLength, msgBuffer, & - msgBufferReduced ) - ! unpack the buffer - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - minLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) - enddo - enddo - enddo - !$omp end do - - ! pack the buffer for maxval - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = maxLayer(iField,iLevel,iRegion) - enddo - enddo - enddo - !$omp end do - call mpas_dmpar_max_real_array(dminfo, bufferLength, msgBuffer, & - msgBufferReduced ) - ! unpack the buffer - !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) - do iRegion=1,nRegions - do iLevel=1,nVertLevels - do iField=1,nFields - bufferAddr = (iRegion -1)*nFields*nVertLevels + & - (iLevel -1)*nFields + iField - maxLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) - enddo - enddo - enddo - !$omp end do - - ! deallocate buffer - !$omp single - deallocate(msgBuffer) - deallocate(msgBufferReduced) - !$omp end single - - ! find min/max with region volume - !$omp do collapse(2) private(iRegion, iField) - do iRegion=1,nRegions - do iField=1,nFields - minField(iField,iRegion) = minval(minLayer(iField,:,iRegion)) - maxField(iField,iRegion) = maxval(maxLayer(iField,:,iRegion)) - enddo - enddo - !$omp end do + if(predefRegions == .true.) then + ! allocate buffer for message passing + bufferAddr=0 + bufferLength=nRegionsTmp*nFields*nVertLevels + !$omp single + allocate(msgBuffer (bufferLength)) + allocate(msgBufferReduced(bufferLength)) + !$omp end single + + ! communication for sums + ! pack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = avgLayer(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_sum_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + + ! compute normalized averages + + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + avgLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + + ! compute vertical sum before layer-by-layer normalization + !$omp do private(iLevel,iField,iRegion,normFactor) + do iRegion=1,nRegionsTmp + do iField=1,nFields + avgField(iField,iRegion) = 0.0_RKIND + do iLevel=1,nVertLevels + avgField(iField,iRegion) = avgField(iField,iRegion) & + + avgLayer(iField,iLevel,iRegion) + enddo + enddo + normFactor = 1.0_RKIND / max(avgField(3,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgField(iField,iRegion) = avgField(iField,iRegion)*normFactor + enddo + ! normalize total region volume by total volume cell area + avgField(3,iRegion) = avgField(3,iRegion) & + / max(avgField(2,iRegion),1.0e-8_RKIND) + ! normalize total volume cell area by total number of cells + avgField(2,iRegion) = avgField(2,iRegion) & + / max(avgField(1,iRegion),1.0e-8_RKIND) + enddo + !$omp end do + + ! normalize averages layer-by-layer + !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + ! normalize all field by total volume in each layer + normFactor = 1.0_RKIND/max(avgLayer(3,iLevel,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgLayer(iField,iLevel,iRegion) = & + avgLayer(iField,iLevel,iRegion) * normFactor + enddo + ! normalize total layer volume by layer area + avgLayer(3,iLevel,iRegion) = avgLayer(3,iLevel,iRegion) & + / max(avgLayer(2,iLevel,iRegion),1.0e-8_RKIND) + ! normalize total layer area by number of cells in region + avgLayer(2,iLevel,iRegion) = avgLayer(2,iLevel,iRegion) & + / max(avgLayer(1,iLevel,iRegion),1.0e-8_RKIND) + enddo + enddo + !$omp end do + + ! now compute min/max for each field, region + ! find global min/max + + ! pack the buffer for minval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = minLayer(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + ! communicate for global min + call mpas_dmpar_min_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + minLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + ! pack the buffer for maxval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = maxLayer(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_max_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegionsTmp + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + maxLayer(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + + ! deallocate buffer + !$omp single + deallocate(msgBuffer) + deallocate(msgBufferReduced) + !$omp end single + endif ! hard-wired regions + + + if(additionalRegion /= '') then + ! allocate buffer for message passing + bufferAddr=0 + bufferLength=nRegions*nFields*nVertLevels + !$omp single + allocate(msgBuffer (bufferLength)) + allocate(msgBufferReduced(bufferLength)) + !$omp end single + + ! communication for sums + ! pack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = avgLayerRF(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_sum_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + + ! compute normalized averages + + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + avgLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + + ! compute vertical sum before layer-by-layer normalization + !$omp do private(iLevel,iField,iRegion,normFactor) + do iRegion=1,nRegions + do iField=1,nFields + avgFieldRF(iField,iRegion) = 0.0_RKIND + do iLevel=1,nVertLevels + avgFieldRF(iField,iRegion) = avgFieldRF(iField,iRegion) & + + avgLayerRF(iField,iLevel,iRegion) + enddo + enddo + normFactor = 1.0_RKIND / max(avgFieldRF(3,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgFieldRF(iField,iRegion) = avgFieldRF(iField,iRegion)*normFactor + enddo + ! normalize total region volume by total volume cell area + avgFieldRF(3,iRegion) = avgFieldRF(3,iRegion) & + / max(avgFieldRF(2,iRegion),1.0e-8_RKIND) + ! normalize total volume cell area by total number of cells + avgFieldRF(2,iRegion) = avgFieldRF(2,iRegion) & + / max(avgFieldRF(1,iRegion),1.0e-8_RKIND) + enddo + !$omp end do + + ! normalize averages layer-by-layer + !$omp do collapse(2) private(iLevel, iRegion, iField, normFactor) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + ! normalize all field by total volume in each layer + normFactor = 1.0_RKIND/max(avgLayerRF(3,iLevel,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgLayerRF(iField,iLevel,iRegion) = & + avgLayerRF(iField,iLevel,iRegion) * normFactor + enddo + ! normalize total layer volume by layer area + avgLayerRF(3,iLevel,iRegion) = avgLayerRF(3,iLevel,iRegion) & + / max(avgLayerRF(2,iLevel,iRegion),1.0e-8_RKIND) + ! normalize total layer area by number of cells in region + avgLayerRF(2,iLevel,iRegion) = avgLayerRF(2,iLevel,iRegion) & + / max(avgLayerRF(1,iLevel,iRegion),1.0e-8_RKIND) + enddo + enddo + !$omp end do + + ! now compute min/max for each field, region + ! find global min/max + + ! pack the buffer for minval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = minLayerRF(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + ! communicate for global min + call mpas_dmpar_min_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + minLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + + ! pack the buffer for maxval + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + msgBuffer(bufferAddr) = maxLayerRF(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_max_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) + ! unpack the buffer + !$omp do collapse(3) private(iRegion,iLevel,iField,bufferAddr) + do iRegion=1,nRegions + do iLevel=1,nVertLevels + do iField=1,nFields + bufferAddr = (iRegion -1)*nFields*nVertLevels + & + (iLevel -1)*nFields + iField + maxLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + enddo + enddo + enddo + !$omp end do + + ! deallocate buffer + !$omp single + deallocate(msgBuffer) + deallocate(msgBufferReduced) + !$omp end single + endif ! region mask files + + if(preDefRegions == .true.) then + ! find min/max with region volume + !$omp do collapse(2) private(iRegion, iField) + do iRegion=1,nRegionsTmp + do iField=1,nFields + minFieldRF(iField,iRegion) = minval(minLayerRF(iField,:,iRegion)) + maxFieldRF(iField,iRegion) = maxval(maxLayerRF(iField,:,iRegion)) + enddo + enddo + !$omp end do + endif ! hard-wired regions + + if(additionalRegion /= '') then + ! find min/max with region volume + !$omp do collapse(2) private(iRegion, iField) + do iRegion=1,nRegions + do iField=1,nFields + minFieldRF(iField,iRegion) = minval(minLayerRF(iField,:,iRegion)) + maxFieldRF(iField,iRegion) = maxval(maxLayerRF(iField,:,iRegion)) + enddo + enddo + !$omp end do + endif ! region mask files + write(stdOutUnit, *) 'maxLayerRF', maxLayerRF + write(stdOutUnit, *) 'maxLayer', maxLayer end subroutine ocn_compute_layer_volume_weighted_averages!}}} !*********************************************************************** @@ -740,9 +1065,9 @@ end subroutine ocn_finalize_layer_volume_weighted_averages!}}} !> \details !> This is a local routine for this analysis member that computes !> a multiplicative mask that is 1 if the point is a member of the -!> desired region and zero otherwise. +!> desired region and zero otherwise. !> NOTE: This routine should be temporary and should be replaced -!> once an MPAS-wide ocean region and mask capability +!> once an MPAS-wide ocean region and mask capability !> are introduced. ! !----------------------------------------------------------------------- @@ -762,7 +1087,7 @@ subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & real(kind=RKIND), dimension(:), intent(in) :: & lonCell, latCell ! lat, lon coordinates for each cell - + !------------------------------------------------------------------ ! ! output variables @@ -779,7 +1104,7 @@ subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & integer :: & iCell, iRegion, &! loop indices - nCells, nRegions, nCellsMax ! array sizes + nCells, nRegionsTmp, nCellsMax ! array sizes real(kind=RKIND) :: & dtr ! degree to radian conversion @@ -787,23 +1112,23 @@ subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & !------------------------------------------------------------------ nCellsMax = size(regionMask, dim=1) - nRegions = size(regionMask, dim=2) + nRegionsTmp = size(regionMask, dim=2) - ! Currently assumes 6 regions, so fail if nRegions not large enough - if (nRegions < 6) then + ! Currently assumes 7 regions, so fail if nRegionsTmp not large enough + if (nRegionsTmp < 7) then !$omp master - write (stderrUnit,*) 'Abort: nRegions inconsistent ' + write (stderrUnit,*) 'Abort: nOceanRegionsTmp inconsistent ' write (stderrUnit,*) ' : in layer volume weighted avg AM' !$omp end master call mpas_dmpar_global_abort(& - 'MPAS-ocean: Abort: nRegions inconsistent in layer vol wgt avg AM') + 'MPAS-ocean: Abort: nOceanRegionsTmp inconsistent in layer vol wgt avg AM') endif ! degrees to radians dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND ! Define masks for each region - !$omp do + !$omp do do iCell=1,nCellsSolve ! Define Arctic region (region 1) @@ -858,6 +1183,9 @@ subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & regionMask(iCell,6,iBlock) = 1.0_RKIND endif + ! global (region 7) + ! do nothing (already initialized to 1) + enddo !$omp end do @@ -865,7 +1193,7 @@ subroutine ocn_LayerVolWgtAvgRegionMask(iBlock, nCellsSolve, & ! extra cells in this block !$omp do collapse(2) - do iRegion=1,nRegions + do iRegion=1,nRegionsTmp do iCell=nCellsSolve+1,nCellsMax regionMask(iCell,iRegion,iBlock) = 0.0_RKIND end do @@ -886,7 +1214,7 @@ end subroutine ocn_LayerVolWgtAvgRegionMask !}}} !> \details !> This is a local routine for this analysis member that computes !> a multiplicative mask that is 1 if a vertical layer is active -!> and is zero otherwise. +!> and is zero otherwise. !> NOTE: This routine should be temporary. An MPAS ocean-wide !> vertical mask will be used eventually. ! @@ -902,7 +1230,7 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & !------------------------------------------------------------------ integer, intent(in) :: & - nCellsSolve ! number of cells owned by this subdomain + nCellsSolve ! number of cells owned by this subdomain integer, intent(in), dimension(:) :: & maxLevelCell ! last active cell in each column @@ -936,13 +1264,13 @@ subroutine ocn_LayerVolWgtAvgVertMask(maxLevelCell, nCellsSolve, & !$omp do private(iLevel,iCell,maxLevel) do iCell=1,nCellsSolve - maxLevel = maxLevelCell(iCell) + maxLevel = maxLevelCell(iCell) do iLevel=1,maxLevel vertMask(iCell,iLevel) = 1.0_RKIND - end do + end do do iLevel=maxLevel+1,nVertLevels vertMask(iCell,iLevel) = 0.0_RKIND - end do + end do end do !$omp end do From 4a5ed7a428131a1093e4aff08b0d6094dc86cb6e Mon Sep 17 00:00:00 2001 From: Anne S Berres Date: Tue, 9 May 2017 17:32:33 -0600 Subject: [PATCH 11/11] Refactored to conform with issue #1305. --- ...egistry_layer_volume_weighted_averages.xml | 204 +++++++++--------- .../mpas_ocn_layer_volume_weighted_averages.F | 112 +++++----- 2 files changed, 158 insertions(+), 158 deletions(-) diff --git a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml index 9fe27ace2a..5aedcd85ef 100644 --- a/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml +++ b/src/core_ocean/analysis_members/Registry_layer_volume_weighted_averages.xml @@ -322,285 +322,285 @@ description="Average relative enstrophy within region volume" /> - - + - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - @@ -623,12 +623,12 @@ - - - - - - + + + + + + diff --git a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F index 48f045c09c..3e9bfb5154 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F +++ b/src/core_ocean/analysis_members/mpas_ocn_layer_volume_weighted_averages.F @@ -235,9 +235,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! pointers to result arrays real (kind=RKIND), dimension(:,:,:), pointer :: & - avgLayer, minLayer, maxLayer, avgLayerRF, minLayerRF, maxLayerRF + avgLayer, minLayer, maxLayer, avgLayerRegionMask, minLayerRegionMask, maxLayerRegionMask real (kind=RKIND), dimension(:,:), pointer :: & - avgField, minField, maxField, avgFieldRF, minFieldRF, maxFieldRF + avgField, minField, maxField, avgFieldRegionMask, minFieldRegionMask, maxFieldRegionMask ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -342,23 +342,23 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ ! get pointers to analysis member arrays to store results call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'minValueWithinOceanLayerRF', & - minLayerRF) + 'minValueWithinOceanLayerRegionMask', & + minLayerRegionMask) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'maxValueWithinOceanLayerRF', & - maxLayerRF) + 'maxValueWithinOceanLayerRegionMask', & + maxLayerRegionMask) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'avgValueWithinOceanLayerRF', & - avgLayerRF) + 'avgValueWithinOceanLayerRegionMask', & + avgLayerRegionMask) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'minValueWithinOceanVolumeRF', & - minFieldRF) + 'minValueWithinOceanVolumeRegionMask', & + minFieldRegionMask) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'maxValueWithinOceanVolumeRF', & - maxFieldRF) + 'maxValueWithinOceanVolumeRegionMask', & + maxFieldRegionMask) call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & - 'avgValueWithinOceanVolumeRF', & - avgFieldRF) + 'avgValueWithinOceanVolumeRegionMask', & + avgFieldRegionMask) !!! region file variables ! if a region file is given and a region is selected, set up region file version @@ -409,9 +409,9 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ avgLayer = 0.0_RKIND minLayer = 1.e20_RKIND maxLayer = -1.e20_RKIND - avgLayerRF = 0.0_RKIND - minLayerRF = 1.e20_RKIND - maxLayerRF = -1.e20_RKIND + avgLayerRegionMask = 0.0_RKIND + minLayerRegionMask = 1.e20_RKIND + maxLayerRegionMask = -1.e20_RKIND ! loop over blocks iBlock = 0 @@ -573,10 +573,10 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iCell=1,nCellsSolve maskTmp = vertMask(iCell,iLevel) * regionCellMasks(iRegion,iCell) if (maskTmp /= 0.0_RKIND) then - minLayerRF(iField,iLevel,iRegion) = min( & - minLayerRF(iField,iLevel,iRegion), fieldTmp(iCell,iField)) - maxLayerRF(iField,iLevel,iRegion) = max( & - maxLayerRF(iField,iLevel,iRegion), fieldTmp(iCell,iField)) + minLayerRegionMask(iField,iLevel,iRegion) = min( & + minLayerRegionMask(iField,iLevel,iRegion), fieldTmp(iCell,iField)) + maxLayerRegionMask(iField,iLevel,iRegion) = max( & + maxLayerRegionMask(iField,iLevel,iRegion), fieldTmp(iCell,iField)) endif end do end do @@ -595,29 +595,29 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iCell=1,nCellsSolve fieldTmp(iCell,1) = vertMask(iCell,iLevel)* & regionCellMasks(iRegion,iCell) - avgLayerRF(1,iLevel,iRegion) = & - avgLayerRF(1,iLevel,iRegion) + fieldTmp(iCell,1) + avgLayerRegionMask(1,iLevel,iRegion) = & + avgLayerRegionMask(1,iLevel,iRegion) + fieldTmp(iCell,1) end do ! cell area do iCell=1,nCellsSolve - avgLayerRF(2,iLevel,iRegion) = & - avgLayerRF(2,iLevel,iRegion) + fieldTmp(iCell,2)* & + avgLayerRegionMask(2,iLevel,iRegion) = & + avgLayerRegionMask(2,iLevel,iRegion) + fieldTmp(iCell,2)* & fieldTmp(iCell,1) end do ! for field three, we want cell volume instead of thickness do iCell=1,nCellsSolve - avgLayerRF(3,iLevel,iRegion) = & - avgLayerRF(3,iLevel,iRegion) + fieldTmp(iCell,3)*& + avgLayerRegionMask(3,iLevel,iRegion) = & + avgLayerRegionMask(3,iLevel,iRegion) + fieldTmp(iCell,3)*& fieldTmp(iCell,1) end do ! for all remaining fields, we want a volume weighted avg do iField=4,nFields do iCell=1,nCellsSolve - avgLayerRF(iField,iLevel,iRegion) = & - avgLayerRF(iField,iLevel,iRegion) + & + avgLayerRegionMask(iField,iLevel,iRegion) = & + avgLayerRegionMask(iField,iLevel,iRegion) + & fieldTmp(iCell,iField)* & fieldTmp(iCell,3) * & fieldTmp(iCell,1) @@ -803,7 +803,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = avgLayerRF(iField,iLevel,iRegion) + msgBuffer(bufferAddr) = avgLayerRegionMask(iField,iLevel,iRegion) enddo enddo enddo @@ -820,7 +820,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - avgLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + avgLayerRegionMask(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo @@ -830,22 +830,22 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ !$omp do private(iLevel,iField,iRegion,normFactor) do iRegion=1,nRegions do iField=1,nFields - avgFieldRF(iField,iRegion) = 0.0_RKIND + avgFieldRegionMask(iField,iRegion) = 0.0_RKIND do iLevel=1,nVertLevels - avgFieldRF(iField,iRegion) = avgFieldRF(iField,iRegion) & - + avgLayerRF(iField,iLevel,iRegion) + avgFieldRegionMask(iField,iRegion) = avgFieldRegionMask(iField,iRegion) & + + avgLayerRegionMask(iField,iLevel,iRegion) enddo enddo - normFactor = 1.0_RKIND / max(avgFieldRF(3,iRegion),1.0e-8_RKIND) + normFactor = 1.0_RKIND / max(avgFieldRegionMask(3,iRegion),1.0e-8_RKIND) do iField=4,nFields - avgFieldRF(iField,iRegion) = avgFieldRF(iField,iRegion)*normFactor + avgFieldRegionMask(iField,iRegion) = avgFieldRegionMask(iField,iRegion)*normFactor enddo ! normalize total region volume by total volume cell area - avgFieldRF(3,iRegion) = avgFieldRF(3,iRegion) & - / max(avgFieldRF(2,iRegion),1.0e-8_RKIND) + avgFieldRegionMask(3,iRegion) = avgFieldRegionMask(3,iRegion) & + / max(avgFieldRegionMask(2,iRegion),1.0e-8_RKIND) ! normalize total volume cell area by total number of cells - avgFieldRF(2,iRegion) = avgFieldRF(2,iRegion) & - / max(avgFieldRF(1,iRegion),1.0e-8_RKIND) + avgFieldRegionMask(2,iRegion) = avgFieldRegionMask(2,iRegion) & + / max(avgFieldRegionMask(1,iRegion),1.0e-8_RKIND) enddo !$omp end do @@ -854,17 +854,17 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iRegion=1,nRegions do iLevel=1,nVertLevels ! normalize all field by total volume in each layer - normFactor = 1.0_RKIND/max(avgLayerRF(3,iLevel,iRegion),1.0e-8_RKIND) + normFactor = 1.0_RKIND/max(avgLayerRegionMask(3,iLevel,iRegion),1.0e-8_RKIND) do iField=4,nFields - avgLayerRF(iField,iLevel,iRegion) = & - avgLayerRF(iField,iLevel,iRegion) * normFactor + avgLayerRegionMask(iField,iLevel,iRegion) = & + avgLayerRegionMask(iField,iLevel,iRegion) * normFactor enddo ! normalize total layer volume by layer area - avgLayerRF(3,iLevel,iRegion) = avgLayerRF(3,iLevel,iRegion) & - / max(avgLayerRF(2,iLevel,iRegion),1.0e-8_RKIND) + avgLayerRegionMask(3,iLevel,iRegion) = avgLayerRegionMask(3,iLevel,iRegion) & + / max(avgLayerRegionMask(2,iLevel,iRegion),1.0e-8_RKIND) ! normalize total layer area by number of cells in region - avgLayerRF(2,iLevel,iRegion) = avgLayerRF(2,iLevel,iRegion) & - / max(avgLayerRF(1,iLevel,iRegion),1.0e-8_RKIND) + avgLayerRegionMask(2,iLevel,iRegion) = avgLayerRegionMask(2,iLevel,iRegion) & + / max(avgLayerRegionMask(1,iLevel,iRegion),1.0e-8_RKIND) enddo enddo !$omp end do @@ -879,7 +879,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = minLayerRF(iField,iLevel,iRegion) + msgBuffer(bufferAddr) = minLayerRegionMask(iField,iLevel,iRegion) enddo enddo enddo @@ -894,7 +894,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - minLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + minLayerRegionMask(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo @@ -907,7 +907,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - msgBuffer(bufferAddr) = maxLayerRF(iField,iLevel,iRegion) + msgBuffer(bufferAddr) = maxLayerRegionMask(iField,iLevel,iRegion) enddo enddo enddo @@ -921,7 +921,7 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ do iField=1,nFields bufferAddr = (iRegion -1)*nFields*nVertLevels + & (iLevel -1)*nFields + iField - maxLayerRF(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) + maxLayerRegionMask(iField,iLevel,iRegion)=msgBufferReduced(bufferAddr) enddo enddo enddo @@ -939,8 +939,8 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ !$omp do collapse(2) private(iRegion, iField) do iRegion=1,nRegionsTmp do iField=1,nFields - minFieldRF(iField,iRegion) = minval(minLayerRF(iField,:,iRegion)) - maxFieldRF(iField,iRegion) = maxval(maxLayerRF(iField,:,iRegion)) + minFieldRegionMask(iField,iRegion) = minval(minLayerRegionMask(iField,:,iRegion)) + maxFieldRegionMask(iField,iRegion) = maxval(maxLayerRegionMask(iField,:,iRegion)) enddo enddo !$omp end do @@ -951,13 +951,13 @@ subroutine ocn_compute_layer_volume_weighted_averages(domain, timeLevel, err)!{{ !$omp do collapse(2) private(iRegion, iField) do iRegion=1,nRegions do iField=1,nFields - minFieldRF(iField,iRegion) = minval(minLayerRF(iField,:,iRegion)) - maxFieldRF(iField,iRegion) = maxval(maxLayerRF(iField,:,iRegion)) + minFieldRegionMask(iField,iRegion) = minval(minLayerRegionMask(iField,:,iRegion)) + maxFieldRegionMask(iField,iRegion) = maxval(maxLayerRegionMask(iField,:,iRegion)) enddo enddo !$omp end do endif ! region mask files - write(stdOutUnit, *) 'maxLayerRF', maxLayerRF + write(stdOutUnit, *) 'maxLayerRegionMask', maxLayerRegionMask write(stdOutUnit, *) 'maxLayer', maxLayer end subroutine ocn_compute_layer_volume_weighted_averages!}}}