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..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 @@ -1,5 +1,5 @@ - + + @@ -314,43 +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 1099e7081a..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 @@ -55,6 +55,12 @@ module ocn_layer_volume_weighted_averages ! !-------------------------------------------------------------------- + ! create a multi-block region mask temporarily + ! this should replaced later by a MPAS ocean-wide region definition + + real (kind=RKIND), dimension(:,:,:), allocatable :: & + regionMask + !*********************************************************************** contains @@ -102,8 +108,72 @@ subroutine ocn_init_layer_volume_weighted_averages(domain, err)!{{{ ! !----------------------------------------------------------------- + type (dm_info) :: dminfo + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: meshPool + + integer, pointer :: & + nCells, &! number of cells in a given block + nCellsSolve, &! number of cells owned by a given block + nRegionsTmp ! number of ocean regions defined + + real (kind=RKIND), dimension(:), pointer :: & + lonCell, latCell ! latitude, longitude of each cell + + integer :: & + iBlock, &! block counter + numBlocks, &! track number of blocks and + maxCells ! largest cell count to size array + + !----------------------------------------------------------------- + err = 0 + ! set highest level pointer + dminfo = domain % dminfo + + ! retrieve number of regions + call mpas_pool_get_dimension(domain % blocklist % dimensions, & + 'nOceanRegionsTmp', nRegionsTmp) + + ! count up the number of blocks and the max size of each domain + numBlocks = 0 + maxCells = 0 + block => 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,nRegionsTmp,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!}}} !*********************************************************************** @@ -111,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. @@ -152,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 @@ -160,13 +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 - 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, avgLayerRegionMask, minLayerRegionMask, maxLayerRegionMask + real (kind=RKIND), dimension(:,:), pointer :: & + avgField, minField, maxField, avgFieldRegionMask, minFieldRegionMask, maxFieldRegionMask ! pointers to data in pools to be analyzed real (kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -182,26 +253,39 @@ 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, nRegionsTmp integer, pointer :: index_temperature, index_salinity integer, dimension(:), pointer :: maxLevelCell - 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 + real (kind=RKIND), dimension(:), pointer :: areaCell + logical, pointer :: predefRegions ! 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 + + !!! 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 ! assume no error err = 0 @@ -209,331 +293,672 @@ 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) + ! 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) + + + call mpas_pool_get_subpool(domain % blocklist % structs, & + 'layerVolumeWeightedAverageAM', & + layerVolumeWeightedAverageAMPool) + + ! 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, & + 'minValueWithinOceanLayerRegionMask', & + minLayerRegionMask) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanLayerRegionMask', & + maxLayerRegionMask) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanLayerRegionMask', & + avgLayerRegionMask) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'minValueWithinOceanVolumeRegionMask', & + minFieldRegionMask) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'maxValueWithinOceanVolumeRegionMask', & + maxFieldRegionMask) + call mpas_pool_get_array(layerVolumeWeightedAverageAMPool, & + 'avgValueWithinOceanVolumeRegionMask', & + avgFieldRegionMask) + + !!! 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 + avgLayerRegionMask = 0.0_RKIND + minLayerRegionMask = 1.e20_RKIND + maxLayerRegionMask = -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, '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 + 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 all ocean regions - do iRegion=1,nOceanRegionsTmp - - ! 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 + + 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 + + if(preDefRegions == .true.) then + ! compute min and max + do iRegion=1,nRegionsTmp + 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,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 + maskTmp = vertMask(iCell,iLevel) * regionCellMasks(iRegion,iCell) + if (maskTmp /= 0.0_RKIND) then + 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 + 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)* & + regionCellMasks(iRegion,iCell) + avgLayerRegionMask(1,iLevel,iRegion) = & + avgLayerRegionMask(1,iLevel,iRegion) + fieldTmp(iCell,1) + end do + + ! cell area + do iCell=1,nCellsSolve + 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 + 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 + avgLayerRegionMask(iField,iLevel,iRegion) = & + avgLayerRegionMask(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 + !$omp single + deallocate(fieldTmp) + deallocate(vertMask) + !$omp end single + + block => block % next + 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 - ! 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,:) - 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, nCellsSolve, workArray, workMask, 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)) + ! 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 - enddo ! iLevel + ! now compute min/max for each field, region + ! find global min/max - end do ! iRegion - kBuffer = 0 + ! 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 - ! 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.) + ! 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) = avgLayerRegionMask(iField,iLevel,iRegion) + enddo + enddo + enddo + !$omp end do + call mpas_dmpar_sum_real_array(dminfo, bufferLength, msgBuffer, & + msgBufferReduced ) - block => block % next - end do + ! compute normalized averages - ! 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 ) - - ! unpack the buffer into intent(out) of this analysis member - kBuffer=0 - do iRegion=1,nOceanRegionsTmp - do iLevel=1,nVertLevels - 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) + ! 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 + avgLayerRegionMask(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 + avgFieldRegionMask(iField,iRegion) = 0.0_RKIND + do iLevel=1,nVertLevels + avgFieldRegionMask(iField,iRegion) = avgFieldRegionMask(iField,iRegion) & + + avgLayerRegionMask(iField,iLevel,iRegion) + enddo enddo - enddo - enddo - - ! 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 - do iLevel=1,nVertLevels - avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField, iRegion) & - + avgValueWithinOceanLayerRegion(iDataField,iLevel,iRegion) - enddo - enddo - do iDataField=4,nDefinedDataFields - avgValueWithinOceanVolumeRegion(iDataField, iRegion) = avgValueWithinOceanVolumeRegion(iDataField, iRegion) & - / max(avgValueWithinOceanVolumeRegion(3,iRegion),1.0e-8_RKIND) - enddo - ! normalize total region volume by total volume cell area - avgValueWithinOceanVolumeRegion(3,iRegion) = avgValueWithinOceanVolumeRegion(3,iRegion) & - / max(avgValueWithinOceanVolumeRegion(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) - enddo - - ! 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 - 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) - enddo - enddo - - ! deallocate buffers - deallocate(workBufferSumReduced) - deallocate(workBufferMinReduced) - deallocate(workBufferMaxReduced) - - contains - - subroutine compute_mask(iLevel, maxLevelCell, nCells, nCellsSolve, iRegion, 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), dimension(:) :: maxLevelCell - real(kind=RKIND), dimension(:), intent(in) :: lonCell, latCell - real(kind=RKIND), dimension(:), intent(out) :: workMask - integer :: iCell - real(kind=RKIND) :: dtr - - dtr = 4.0_RKIND*atan(1.0_RKIND) / 180.0_RKIND - workMask(:) = 0.0_RKIND - do iCell=1,nCellsSolve - if(iLevel.le.maxLevelCell(iCell)) workMask(iCell) = 1.0_RKIND - 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 - - end subroutine compute_mask - - - subroutine compute_statistics(nDefinedDataFields, 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 - 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 - 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 + normFactor = 1.0_RKIND / max(avgFieldRegionMask(3,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgFieldRegionMask(iField,iRegion) = avgFieldRegionMask(iField,iRegion)*normFactor + enddo + ! normalize total region volume by total volume cell area + avgFieldRegionMask(3,iRegion) = avgFieldRegionMask(3,iRegion) & + / max(avgFieldRegionMask(2,iRegion),1.0e-8_RKIND) + ! normalize total volume cell area by total number of cells + avgFieldRegionMask(2,iRegion) = avgFieldRegionMask(2,iRegion) & + / max(avgFieldRegionMask(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(avgLayerRegionMask(3,iLevel,iRegion),1.0e-8_RKIND) + do iField=4,nFields + avgLayerRegionMask(iField,iLevel,iRegion) = & + avgLayerRegionMask(iField,iLevel,iRegion) * normFactor + enddo + ! normalize total layer volume by layer area + 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 + avgLayerRegionMask(2,iLevel,iRegion) = avgLayerRegionMask(2,iLevel,iRegion) & + / max(avgLayerRegionMask(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) = minLayerRegionMask(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 + minLayerRegionMask(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) = maxLayerRegionMask(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 + maxLayerRegionMask(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 + minFieldRegionMask(iField,iRegion) = minval(minLayerRegionMask(iField,:,iRegion)) + maxFieldRegionMask(iField,iRegion) = maxval(maxLayerRegionMask(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 + 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, *) 'maxLayerRegionMask', maxLayerRegionMask + write(stdOutUnit, *) 'maxLayer', maxLayer end subroutine ocn_compute_layer_volume_weighted_averages!}}} !*********************************************************************** @@ -630,6 +1055,229 @@ 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, nRegionsTmp, nCellsMax ! array sizes + + real(kind=RKIND) :: & + dtr ! degree to radian conversion + + !------------------------------------------------------------------ + + nCellsMax = size(regionMask, dim=1) + nRegionsTmp = size(regionMask, dim=2) + + ! Currently assumes 7 regions, so fail if nRegionsTmp not large enough + if (nRegionsTmp < 7) then + !$omp master + write (stderrUnit,*) 'Abort: nOceanRegionsTmp inconsistent ' + write (stderrUnit,*) ' : in layer volume weighted avg AM' + !$omp end master + call mpas_dmpar_global_abort(& + '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 + 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 + + ! global (region 7) + ! do nothing (already initialized to 1) + + 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,nRegionsTmp + 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 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 @@ + + - + + + + + + + + +