From 942333cea6b28cf00e3981467594a4771b9900a1 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:30:50 -0600 Subject: [PATCH 01/46] minLevel* changes to ocn_tracer_exponential_decay --- src/core_ocean/shared/mpas_ocn_tendency.F | 3 ++- src/core_ocean/shared/mpas_ocn_tracer_exponential_decay.F | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index dde9ac3d3c..4e1b281d80 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -571,6 +571,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing call mpas_pool_get_array(forcingPool, 'penetrativeTemperatureFlux', penetrativeTemperatureFlux) call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed) call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'latCell', latCell) @@ -785,7 +786,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing modifiedGroupName = trim(groupItr % memberName) // "ExponentialDecayRate" call mpas_pool_get_array(tracersExponentialDecayFieldsPool, trim(modifiedGroupName), & tracerGroupExponentialDecayRate) - call ocn_tracer_exponential_decay_compute(nTracerGroup, nCells, maxLevelCell, layerThickness, & + call ocn_tracer_exponential_decay_compute(nTracerGroup, nCells, minLevelCell, maxLevelCell, layerThickness, & tracerGroup, tracerGroupExponentialDecayRate, tracerGroupTend, err) call mpas_timer_stop("exponential decay " // trim(groupItr % memberName)) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_exponential_decay.F b/src/core_ocean/shared/mpas_ocn_tracer_exponential_decay.F index 63379c8b2a..676f702571 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_exponential_decay.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_exponential_decay.F @@ -65,8 +65,8 @@ module ocn_tracer_exponential_decay ! !----------------------------------------------------------------------- - subroutine ocn_tracer_exponential_decay_compute(nTracers, nCellsSolve, maxLevelCell, layerThickness, tracers, & !{{{ - tracersExponentialDecayRate, tracer_tend, err) + subroutine ocn_tracer_exponential_decay_compute(nTracers, nCellsSolve, minLevelCell, maxLevelCell, & !{{{ + layerThickness, tracers, tracersExponentialDecayRate, tracer_tend, err) !----------------------------------------------------------------- ! @@ -76,7 +76,7 @@ subroutine ocn_tracer_exponential_decay_compute(nTracers, nCellsSolve, maxLevelC ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell real (kind=RKIND), dimension(:), intent(in) :: & tracersExponentialDecayRate @@ -122,7 +122,7 @@ subroutine ocn_tracer_exponential_decay_compute(nTracers, nCellsSolve, maxLevelC !$omp parallel !$omp do schedule(runtime) private(iLevel, iTracer) do iCell=1,nCellsSolve - do iLevel=1,maxLevelCell(iCell) + do iLevel=minLevelCell(iCell),maxLevelCell(iCell) do iTracer=1,nTracers tracer_tend(iTracer,iLevel,iCell) = tracer_tend(iTracer,iLevel,iCell) & - ( layerThickness(iLevel,iCell) & From 531896a163094a8d5837d19ea1c67f6d1727fdaa Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:38:37 -0600 Subject: [PATCH 02/46] minLevel* changes to ocn_tracer_DMS --- src/core_ocean/shared/mpas_ocn_tendency.F | 4 ++-- src/core_ocean/shared/mpas_ocn_tracer_DMS.F | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 4e1b281d80..4961e0a7bc 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -702,11 +702,11 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel) nTracersEcosys = size(ecosysTracers, dim=1) call ocn_tracer_DMS_compute(activeTracers, tracerGroup, nTracerGroup, ecosysTracers, & - nTracersEcosys, forcingPool, nCellsSolve, maxLevelCell, & + nTracersEcosys, forcingPool, nCellsSolve, minLevelCell, maxLevelCell, & nVertLevels, layerThickness, indexTemperature, indexSalinity, tracerGroupTend, err) call ocn_tracer_DMS_surface_flux_compute(activeTracers, tracerGroup, forcingPool, & - nTracerGroup, nCellsSolve, zMid, indexTemperature, indexSalinity, tracerGroupSurfaceFlux, & + nTracerGroup, nCellsSolve, zMid, minLevelCell, indexTemperature, indexSalinity, tracerGroupSurfaceFlux, & tracerGroupSurfaceFluxRemoved, err) endif diff --git a/src/core_ocean/shared/mpas_ocn_tracer_DMS.F b/src/core_ocean/shared/mpas_ocn_tracer_DMS.F index 5d91827ae5..0076c867fb 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_DMS.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_DMS.F @@ -97,7 +97,7 @@ module ocn_tracer_DMS !----------------------------------------------------------------------- subroutine ocn_tracer_DMS_compute(activeTracers, DMSTracers, nTracersDMS, ecosysTracers, nTracersEcosys, & - forcingPool, nCellsSolve, maxLevelCell, & + forcingPool, nCellsSolve, minLevelCell, maxLevelCell, & nVertLevels, layerThickness, indexTemperature, indexSalinity, DMSTracersTend, err)!{{{ !----------------------------------------------------------------- @@ -108,7 +108,7 @@ subroutine ocn_tracer_DMS_compute(activeTracers, DMSTracers, nTracersDMS, ecosys ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell ! two dimensional arrays real (kind=RKIND), dimension(:,:), intent(in) :: & @@ -172,18 +172,18 @@ subroutine ocn_tracer_DMS_compute(activeTracers, DMSTracers, nTracersDMS, ecosys numColumns = 1 column = 1 - iLevelSurface = 1 + iLevelSurface = minLevelCell(iCell) !DWJ 08/05/2016: This loop needs OpenMP added to it. do iCell=1,nCellsSolve - DMS_input%number_of_active_levels(column) = maxLevelCell(iCell) + DMS_input%number_of_active_levels(column) = maxLevelCell(iCell) - minLevelCell(iCell) + 1 DMS_forcing%ShortWaveFlux_surface(column) = shortWaveHeatFlux(iCell) DMS_forcing%SST(column) = activeTracers(indexTemperature,iLevelSurface,iCell) DMS_forcing%SSS(column) = activeTracers(indexSalinity,iLevelSurface,iCell) - do iLevel=1,maxLevelCell(iCell) - DMS_input%cell_thickness(iLevel,column) = layerThickness(iLevel,iCell)*convertLengthScale + do iLevel=1,DMS_input%number_of_active_levels(column) + DMS_input%cell_thickness(iLevel,column) = layerThickness(iLevel + minLevelCell(iCell) - 1,iCell)*convertLengthScale DMS_input%DMS_tracers(iLevel,column,DMS_indices%dms_ind) = DMSTracers(dmsIndices%dms_ind,iLevel,iCell) DMS_input%DMS_tracers(iLevel,column,DMS_indices%dmsp_ind) = DMSTracers(dmsIndices%dmsp_ind,iLevel,iCell) @@ -207,7 +207,7 @@ subroutine ocn_tracer_DMS_compute(activeTracers, DMSTracers, nTracersDMS, ecosys DMS_output, DMS_diagnostic_fields, nVertLevels, & numColumnsMax, numColumns) - do iLevel=1,maxLevelCell(iCell) + do iLevel=1,DMS_input%number_of_active_levels(column) DMSTracersTend(dmsIndices%dms_ind,iLevel,iCell) = DMSTracersTend(dmsIndices%dms_ind,iLevel,iCell) & + DMS_output%DMS_tendencies(iLevel,column,DMS_indices%dms_ind)*layerThickness(iLevel,iCell) @@ -237,7 +237,7 @@ end subroutine ocn_tracer_DMS_compute!}}} !----------------------------------------------------------------------- subroutine ocn_tracer_DMS_surface_flux_compute(activeTracers, DMSTracers, forcingPool, & - nTracers, nCellsSolve, zMid, indexTemperature, indexSalinity, DMSSurfaceFlux, DMSSurfaceFluxRemoved, err)!{{{ + nTracers, nCellsSolve, zMid, minLevelCell, indexTemperature, indexSalinity, DMSSurfaceFlux, DMSSurfaceFluxRemoved, err)!{{{ !----------------------------------------------------------------- ! @@ -260,6 +260,7 @@ subroutine ocn_tracer_DMS_surface_flux_compute(activeTracers, DMSTracers, forcin ! scalars integer, intent(in) :: nTracers, nCellsSolve, indexTemperature, indexSalinity + integer, dimension(:), intent(in) :: minLevelCell type (mpas_pool_type), intent(inout) :: forcingPool @@ -349,11 +350,11 @@ subroutine ocn_tracer_DMS_surface_flux_compute(activeTracers, DMSTracers, forcin numColumns = 1 column = 1 - iLevelSurface = 1 !DWJ 08/05/2016: This loop needs OpenMP added to it. do iCell=1,nCellsSolve + iLevelSurface = minLevelCell(iCell) DMS_forcing%surfacePressure(column) = atmosphericPressure(iCell)*PascalsToAtmospheres DMS_forcing%iceFraction(column) = iceFraction(iCell) !maltrud assume for now that if there is any land ice, it is all land ice From 2e698a37deb4010007f89fdb93bc34dea8bfb30e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:49:14 -0600 Subject: [PATCH 03/46] minLevel* changes to ocn_tracer_ecosys --- src/core_ocean/shared/mpas_ocn_tendency.F | 4 ++-- src/core_ocean/shared/mpas_ocn_tracer_ecosys.F | 17 ++++++++++------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 4961e0a7bc..33c4051fc8 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -687,11 +687,11 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing if ( trim(groupItr % memberName) == 'ecosysTracers' ) then call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, timeLevel) call ocn_tracer_ecosys_compute(activeTracers, tracerGroup, forcingPool, nTracerGroup, & - nCellsSolve, latCell, maxLevelCell, nVertLevels, layerThickness, zMid, indexTemperature, & + nCellsSolve, latCell, minLevelCell, maxLevelCell, nVertLevels, layerThickness, zMid, indexTemperature, & indexSalinity, tracerGroupTend, err) call ocn_tracer_ecosys_surface_flux_compute(activeTracers, tracerGroup, forcingPool, & - nTracerGroup, nCellsSolve, zMid, indexTemperature, indexSalinity, tracerGroupSurfaceFlux, err) + nTracerGroup, nCellsSolve, zMid, minLevelCell, indexTemperature, indexSalinity, tracerGroupSurfaceFlux, err) endif ! diff --git a/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F b/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F index 3cb53fca47..9c2c166077 100755 --- a/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F @@ -101,7 +101,7 @@ module ocn_tracer_ecosys !----------------------------------------------------------------------- subroutine ocn_tracer_ecosys_compute(activeTracers, ecosysTracers, forcingPool, nTracers, nCellsSolve, & - latCell, maxLevelCell, nVertLevels, layerThickness, zMid, indexTemperature, indexSalinity, ecosysTracersTend, err)!{{{ + latCell, minLevelCell, maxLevelCell, nVertLevels, layerThickness, zMid, indexTemperature, indexSalinity, ecosysTracersTend, err)!{{{ !----------------------------------------------------------------- ! @@ -111,7 +111,7 @@ subroutine ocn_tracer_ecosys_compute(activeTracers, ecosysTracers, forcingPool, ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell real (kind=RKIND), dimension(:), intent(in) :: & latCell @@ -491,8 +491,9 @@ subroutine ocn_tracer_ecosys_compute(activeTracers, ecosysTracers, forcingPool, ! convert dust flux from kg/m2/s to g/cm2/s BGC_forcing%dust_FLUX_IN(column) = dust_FLUX_IN(iCell)*0.1_RKIND BGC_forcing%ShortWaveFlux_surface(column) = shortWaveHeatFlux(iCell) - zTop = 0.0_RKIND - do iLevel=1,maxLevelCell(iCell) + zTop = 0.0_RKIND ! doesn't appear to be used but may need to changed if ecosys below ice shelves + ! iLevel may have to be adjusted in relation to minLevelCell on the LHS + do iLevel=minLevelCell(iCell),maxLevelCell(iCell) BGC_input%PotentialTemperature(iLevel,column) = activeTracers(indexTemperature,iLevel,iCell) BGC_input%Salinity(iLevel,column) = activeTracers(indexSalinity,iLevel,iCell) BGC_input%cell_center_depth(iLevel,column) = -1.0_RKIND*zMid(iLevel,iCell)*convertLengthMKStoCGS @@ -577,7 +578,7 @@ subroutine ocn_tracer_ecosys_compute(activeTracers, ecosysTracers, forcingPool, BGC_output, BGC_diagnostic_fields, nVertLevels, & numColumnsMax, numColumns, config_ecosys_atm_alt_co2_use_eco) - do iLevel=1,maxLevelCell(iCell) + do iLevel=minLevelCell(iCell),maxLevelCell(iCell) ecosysTracersTend(ecosysIndices%po4_ind,iLevel,iCell) = & ecosysTracersTend(ecosysIndices%po4_ind,iLevel,iCell) & @@ -966,7 +967,7 @@ end subroutine ocn_tracer_ecosys_compute!}}} !----------------------------------------------------------------------- subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, forcingPool, & - nTracers, nCellsSolve, zMid, indexTemperature, indexSalinity, ecosysSurfaceFlux, err)!{{{ + nTracers, nCellsSolve, zMid, minLevelCell, indexTemperature, indexSalinity, ecosysSurfaceFlux, err)!{{{ !----------------------------------------------------------------- ! @@ -988,6 +989,7 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, ! scalars integer, intent(in) :: nTracers, nCellsSolve, indexTemperature, indexSalinity + integer, dimension(:), intent(in) :: minLevelCell type (mpas_pool_type), intent(inout) :: forcingPool @@ -1168,10 +1170,10 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, numColumns = 1 column = 1 - iLevelSurface = 1 !DWJ 08/05/2016: This loop needs OpenMP added to it. do iCell=1,nCellsSolve + iLevelSurface = minLevelCell(iCell) ! BGC_forcing%surfacePressure(column) = atmosphericPressure(iCell)*PascalsToAtmospheres BGC_forcing%surfacePressure(column) = 1.0_RKIND BGC_forcing%iceFraction(column) = iceFraction(iCell) @@ -1188,6 +1190,7 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, BGC_forcing%SST(column) = activeTracers(indexTemperature,iLevelSurface,iCell) BGC_forcing%SSS(column) = activeTracers(indexSalinity,iLevelSurface,iCell) + ! the 1 indices may need to be changed to iLevelSurface BGC_input%BGC_tracers(1,column,BGC_indices%dic_ind) = ecosysTracers(ecosysIndices%dic_ind,1,iCell) BGC_input%BGC_tracers(1,column,BGC_indices%dic_alt_co2_ind) = ecosysTracers(ecosysIndices%dic_alt_co2_ind,1,iCell) BGC_input%BGC_tracers(1,column,BGC_indices%alk_ind) = ecosysTracers(ecosysIndices%alk_ind,1,iCell) From 11c212158bb1cf1fc364d26ff5f891b97ce2022c Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:10:52 -0600 Subject: [PATCH 04/46] minLevel* changes to ocn_tracer_MacroMolecules --- src/core_ocean/shared/mpas_ocn_tendency.F | 2 +- .../shared/mpas_ocn_tracer_MacroMolecules.F | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 33c4051fc8..05b5865f30 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -718,7 +718,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing call mpas_pool_get_array(tracersPool, 'ecosysTracers', ecosysTracers, timeLevel) nTracersEcosys = size(ecosysTracers, dim=1) call ocn_tracer_MacroMolecules_compute(tracerGroup, nTracerGroup, ecosysTracers, nTracersEcosys, forcingPool, & - nCellsSolve, maxLevelCell, nVertLevels, layerThickness, & + nCellsSolve, minLevelCell, maxLevelCell, nVertLevels, layerThickness, & tracerGroupTend, err) call ocn_tracer_MacroMolecules_surface_flux_compute(activeTracers, tracerGroup, forcingPool, & diff --git a/src/core_ocean/shared/mpas_ocn_tracer_MacroMolecules.F b/src/core_ocean/shared/mpas_ocn_tracer_MacroMolecules.F index 8597c07258..76f8c4906f 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_MacroMolecules.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_MacroMolecules.F @@ -91,7 +91,7 @@ module ocn_tracer_MacroMolecules subroutine ocn_tracer_MacroMolecules_compute(MacroMoleculesTracers, nTracersMacroMolecules, & ecosysTracers, nTracersEcosys, forcingPool, & - nCellsSolve, maxLevelCell, nVertLevels, layerThickness, MacroMoleculesTracersTend, err)!{{{ + nCellsSolve, minLevelCell, maxLevelCell, nVertLevels, layerThickness, MacroMoleculesTracersTend, err)!{{{ !----------------------------------------------------------------- ! @@ -101,7 +101,7 @@ subroutine ocn_tracer_MacroMolecules_compute(MacroMoleculesTracers, nTracersMacr ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell ! two dimensional arrays real (kind=RKIND), dimension(:,:), intent(in) :: & @@ -158,9 +158,9 @@ subroutine ocn_tracer_MacroMolecules_compute(MacroMoleculesTracers, nTracersMacr column = 1 !DWJ 08/05/2016: This loop needs OpenMP added to it do iCell=1,nCellsSolve - MacroMolecules_input%number_of_active_levels(column) = maxLevelCell(iCell) - do iLevel=1,maxLevelCell(iCell) - MacroMolecules_input%cell_thickness(iLevel,column) = layerThickness(iLevel,iCell)*convertLengthScale + MacroMolecules_input%number_of_active_levels(column) = maxLevelCell(iCell) - minLevelCell(iCell) + 1 + do iLevel=1,MacroMolecules_input%number_of_active_levels(column) + MacroMolecules_input%cell_thickness(iLevel,column) = layerThickness(iLevel + minLevelCell(iCell) - 1,iCell)*convertLengthScale MacroMolecules_input%MACROS_tracers(iLevel,column,MacroMolecules_indices%prot_ind) = & MacroMoleculesTracers(macrosIndices%prot_ind,iLevel,iCell) @@ -186,7 +186,7 @@ subroutine ocn_tracer_MacroMolecules_compute(MacroMoleculesTracers, nTracersMacr MacroMolecules_output, MacroMolecules_diagnostic_fields, nVertLevels, & numColumnsMax, numColumns) - do iLevel=1,maxLevelCell(iCell) + do iLevel=1,MacroMolecules_input%number_of_active_levels(column) MacroMoleculesTracersTend(macrosIndices%prot_ind,iLevel,iCell) = & MacroMoleculesTracersTend(macrosIndices%prot_ind,iLevel,iCell) & From 517799574219886a27b93871c2e7f47bc60a9613 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 16:03:49 -0600 Subject: [PATCH 05/46] minLevel* changes to ocn_tracer_ideal_age --- src/core_ocean/shared/mpas_ocn_tendency.F | 2 +- src/core_ocean/shared/mpas_ocn_tracer_ideal_age.F | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 05b5865f30..f585e16c63 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -804,7 +804,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing call mpas_pool_get_subpool(forcingPool, 'tracersIdealAgeFields', tracersIdealAgeFieldsPool) modifiedGroupName = trim(groupItr % memberName) // "IdealAgeMask" call mpas_pool_get_array(tracersIdealAgeFieldsPool, trim(modifiedGroupName), tracerGroupIdealAgeMask) - call ocn_tracer_ideal_age_compute(nTracerGroup, nCells, maxLevelCell, layerThickness, & + call ocn_tracer_ideal_age_compute(nTracerGroup, nCells, minLevelCell, maxLevelCell, layerThickness, & tracerGroupIdealAgeMask, tracerGroup, tracerGroupTend, err) call mpas_timer_stop("ideal age " // trim(groupItr % memberName)) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_ideal_age.F b/src/core_ocean/shared/mpas_ocn_tracer_ideal_age.F index 8188601981..2b299553bc 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_ideal_age.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_ideal_age.F @@ -65,7 +65,7 @@ module ocn_tracer_ideal_age ! !----------------------------------------------------------------------- - subroutine ocn_tracer_ideal_age_compute(nTracers, nCellsSolve, maxLevelCell, layerThickness, & + subroutine ocn_tracer_ideal_age_compute(nTracers, nCellsSolve, minLevelCell, maxLevelCell, layerThickness, & idealAgeMask, tracers, tracer_tend, err)!{{{ !----------------------------------------------------------------- @@ -76,7 +76,7 @@ subroutine ocn_tracer_ideal_age_compute(nTracers, nCellsSolve, maxLevelCell, lay ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell ! two dimensional arrays real (kind=RKIND), dimension(:,:), intent(in) :: & @@ -121,7 +121,7 @@ subroutine ocn_tracer_ideal_age_compute(nTracers, nCellsSolve, maxLevelCell, lay !$omp parallel !$omp do schedule(runtime) private(iLevel, iTracer) do iCell=1,nCellsSolve - do iLevel=1,maxLevelCell(iCell) + do iLevel=minLevelCell(iCell),maxLevelCell(iCell) do iTracer=1,nTracers ! zero tracers at surface to zero where idealAgeMask == zero ! idealAgeMask should be equal to 1.0 elsewhere From 2280a9869b1e1b8364d6cf073c5db3bf66e6ef1b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:23:46 -0600 Subject: [PATCH 06/46] minLevel* changes to ocn_tracer_interior_restoring --- src/core_ocean/shared/mpas_ocn_tendency.F | 2 +- src/core_ocean/shared/mpas_ocn_tracer_interior_restoring.F | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index f585e16c63..0233d755f5 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -769,7 +769,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, meshPool, swForcing modifiedGroupName = trim(groupItr % memberName) // "InteriorRestoringValue" call mpas_pool_get_array(tracersInteriorRestoringFieldsPool, trim(modifiedGroupName), & tracerGroupInteriorRestoringValue) - call ocn_tracer_interior_restoring_compute(nTracerGroup, nCells, maxLevelCell, layerThickness, & + call ocn_tracer_interior_restoring_compute(nTracerGroup, nCells, minLevelCell, maxLevelCell, layerThickness, & tracerGroup, tracerGroupInteriorRestoringRate, tracerGroupInteriorRestoringValue, tracerGroupTend, err) call mpas_timer_stop("interior_restoring_" // trim(groupItr % memberName)) endif diff --git a/src/core_ocean/shared/mpas_ocn_tracer_interior_restoring.F b/src/core_ocean/shared/mpas_ocn_tracer_interior_restoring.F index 5b5c49f2a9..bcd27f580c 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_interior_restoring.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_interior_restoring.F @@ -65,7 +65,7 @@ module ocn_tracer_interior_restoring ! !----------------------------------------------------------------------- - subroutine ocn_tracer_interior_restoring_compute(nTracers, nCellsSolve, maxLevelCell, layerThickness, & + subroutine ocn_tracer_interior_restoring_compute(nTracers, nCellsSolve, minLevelCell, maxLevelCell, layerThickness, & tracers, tracersInteriorRestoringRate, tracersInteriorRestoringValue, tracer_tend, err)!{{{ !----------------------------------------------------------------- @@ -76,7 +76,7 @@ subroutine ocn_tracer_interior_restoring_compute(nTracers, nCellsSolve, maxLevel ! one dimensional arrays integer, dimension(:), intent(in) :: & - maxLevelCell + minLevelCell, maxLevelCell ! two dimensional arrays real (kind=RKIND), dimension(:,:), intent(in) :: & @@ -121,7 +121,7 @@ subroutine ocn_tracer_interior_restoring_compute(nTracers, nCellsSolve, maxLevel !$omp parallel !$omp do schedule(runtime) private(iLevel, iTracer) do iCell=1,nCellsSolve - do iLevel=1,maxLevelCell(iCell) + do iLevel=minLevelCell(iCell),maxLevelCell(iCell) do iTracer=1,nTracers tracer_tend(iTracer, iLevel, iCell) = tracer_tend(iTracer, iLevel, iCell) - layerThickness(iLevel,iCell) & * ( tracers(iTracer, iLevel, iCell) & From 88ae8cca0217fe2d4358bc6f9f11867c6d271abd Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:18:36 -0600 Subject: [PATCH 07/46] minLevel* changes to ocn_filter_btr_mode_tend_vel --- src/core_ocean/shared/mpas_ocn_diagnostics.F | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics.F b/src/core_ocean/shared/mpas_ocn_diagnostics.F index 21055eda0c..06101a3312 100644 --- a/src/core_ocean/shared/mpas_ocn_diagnostics.F +++ b/src/core_ocean/shared/mpas_ocn_diagnostics.F @@ -1855,16 +1855,17 @@ subroutine ocn_filter_btr_mode_tend_vel(tendPool, statePool, meshPool, timeLevel ! thicknessSum is initialized outside the loop because on land boundaries ! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a ! nonzero value to avoid a NaN. - normalThicknessFluxSum = layerThickEdge(1,iEdge) * tend_normalVelocity(1,iEdge) - thicknessSum = layerThickEdge(1,iEdge) + normalThicknessFluxSum = layerThickEdge(maxLevelEdgeBot(iEdge),iEdge) * & + tend_normalVelocity(maxLevelEdgeBot(iEdge),iEdge) + thicknessSum = layerThickEdge(maxLevelEdgeBot(iEdge),iEdge) - do k = 2, maxLevelEdgeTop(iEdge) + do k = maxLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) normalThicknessFluxSum = normalThicknessFluxSum + layerThickEdge(k,iEdge) * tend_normalVelocity(k,iEdge) thicknessSum = thicknessSum + layerThickEdge(k,iEdge) enddo vertSum = normalThicknessFluxSum / thicknessSum - do k = 1, maxLevelEdgeTop(iEdge) + do k = maxLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) tend_normalVelocity(k,iEdge) = tend_normalVelocity(k,iEdge) - vertSum enddo enddo ! iEdge From 2f0c2b177f9b1cd38d2c189aa05cc00596a4d586 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:20:44 -0600 Subject: [PATCH 08/46] minLevel* changes to ocn_tend_freq_filtered_thickness --- src/core_ocean/shared/mpas_ocn_tendency.F | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 0233d755f5..85818bd987 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -1035,7 +1035,9 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, meshPool, timeL call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', maxLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) @@ -1070,15 +1072,15 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, meshPool, timeL do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) flux = layerThickEdge(k, iEdge) * normalVelocity(k, iEdge) * dvEdge(iEdge) * edgeSignOnCell(i, iCell) * invAreaCell div_hu(k) = div_hu(k) - flux div_hu_btr = div_hu_btr - flux end do end do - totalThickness = sum(layerThickness(1:maxLevelCell(iCell),iCell)) - do k = 1, maxLevelCell(iCell) + totalThickness = sum(layerThickness(minLevelCell(iCell):maxLevelCell(iCell),iCell)) + do k = minLevelCell(iCell), maxLevelCell(iCell) tend_lowFreqDivergence(k,iCell) = & -2.0 * pii / thickness_filter_timescale_sec & From 4833ca54e58d69bc5f26e062a150a8aadfca8bd6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:18:10 -0600 Subject: [PATCH 09/46] minLevel* changes to ocn_tracer_short_wave_absorption_jerlov --- .../mpas_ocn_tracer_short_wave_absorption_jerlov.F | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_jerlov.F b/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_jerlov.F index 89ee3baeb3..5954b53e73 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_jerlov.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_jerlov.F @@ -128,7 +128,7 @@ subroutine ocn_tracer_short_wave_absorption_jerlov_tend(meshPool, forcingPool, i integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND) :: depth, fluxRemaining real (kind=RKIND), dimension(:), pointer :: refBottomDepth @@ -139,6 +139,7 @@ subroutine ocn_tracer_short_wave_absorption_jerlov_tend(meshPool, forcingPool, i call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth) @@ -154,9 +155,9 @@ subroutine ocn_tracer_short_wave_absorption_jerlov_tend(meshPool, forcingPool, i !$omp firstprivate(weights) #endif do iCell = 1, nCells - depth = 0.0_RKIND + depth = 0.0_RKIND ! depth relative to minLevelCell(iCell) fluxRemaining = 1.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) depth = depth + layerThickness(k, iCell) call ocn_get_jerlov_fraction(depth, weights(k+1)) @@ -174,13 +175,13 @@ subroutine ocn_tracer_short_wave_absorption_jerlov_tend(meshPool, forcingPool, i end if depth = 0.0_RKIND - do k=1,maxLevelCell(iCell) + do k=minLevelCell(iCell),maxLevelCell(iCell) depth = depth + layerThickness(k,iCell) if(depth > abs(config_surface_buoyancy_depth)) exit enddo - if(k == maxLevelCell(iCell) .or. k == 1) then - depLev=2 + if(k == maxLevelCell(iCell) .or. k == minLevelCell(iCell)) then + depLev=minLevelCell(iCell)+1 else depLev=k endif From 97d835e78294fa21d002a063950d7a6ffae4b3a9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:21:19 -0600 Subject: [PATCH 10/46] minLevel* changes to ocn_tracer_short_wave_absorption_variable --- ...mpas_ocn_tracer_short_wave_absorption_variable.F | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_variable.F b/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_variable.F index bd1343266e..f458141eac 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_variable.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_short_wave_absorption_variable.F @@ -127,7 +127,7 @@ subroutine ocn_tracer_short_wave_absorption_variable_tend(meshPool, swForcingPoo integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND) :: depth real (kind=RKIND), dimension(:), pointer :: refBottomDepth @@ -141,12 +141,13 @@ subroutine ocn_tracer_short_wave_absorption_variable_tend(meshPool, swForcingPoo call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth) allocate(weights(nVertLevels+1)) weights = 0.0_RKIND - weights(1) = 1.0_RKIND + weights(minLevelCell(iCell)) = 1.0_RKIND Avals(:)=0.0_RKIND Kvals(:)=0.0_RKIND @@ -168,7 +169,7 @@ subroutine ocn_tracer_short_wave_absorption_variable_tend(meshPool, swForcingPoo call ocn_get_os00_coeffs(chlorophyllA(iCell),zenithAngle(iCell),cloudRatio,Avals, Kvals) - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) depth = depth + layerThickness(k, iCell) call ocn_get_variable_sw_fraction(depth, weights(k+1), Avals, Kvals) @@ -186,13 +187,13 @@ subroutine ocn_tracer_short_wave_absorption_variable_tend(meshPool, swForcingPoo end if depth = 0.0_RKIND - do k=1,maxLevelCell(iCell) + do k=minLevelCell(iCell),maxLevelCell(iCell) depth = depth + layerThickness(k,iCell) if(depth > abs(config_surface_buoyancy_depth)) exit enddo - if(k == maxLevelCell(iCell) .or. k == 1) then - depLev=2 + if(k == maxLevelCell(iCell) .or. k == minLevelCell(iCell)) then + depLev=minLevelCell(iCell)+1 else depLev=k endif From d0618e64fdbf1d27c90d4b0d0cb937856eb8c2fa Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 11:56:07 -0600 Subject: [PATCH 11/46] minLevel* changes to ocn_frazil_forcing --- .../shared/mpas_ocn_frazil_forcing.F | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_frazil_forcing.F b/src/core_ocean/shared/mpas_ocn_frazil_forcing.F index 74bd39af66..1162434380 100644 --- a/src/core_ocean/shared/mpas_ocn_frazil_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_frazil_forcing.F @@ -167,7 +167,7 @@ subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThickn integer :: iCell, k, nCells integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND), dimension(:,:), pointer :: frazilLayerThicknessTendency err = 0 @@ -177,6 +177,7 @@ subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThickn call mpas_timer_start("frazil thickness tendency") call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(forcingPool, 'frazilLayerThicknessTendency', frazilLayerThicknessTendency) @@ -187,7 +188,7 @@ subroutine ocn_frazil_forcing_layer_thickness(meshPool, forcingPool, layerThickn !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessTend(k,iCell) = layerThicknessTend(k,iCell) + frazilLayerThicknessTendency(k,iCell) end do end do @@ -248,7 +249,7 @@ subroutine ocn_frazil_forcing_active_tracers(meshPool, tracersPool, forcingPool, integer, pointer :: indexTemperature integer, pointer :: indexSalinity integer, pointer, dimension(:) :: nCellsArray - integer, pointer, dimension(:) :: maxLevelCell + integer, pointer, dimension(:) :: minLevelCell, maxLevelCell real (kind=RKIND), dimension(:,:), pointer :: frazilTemperatureTendency real (kind=RKIND), dimension(:,:), pointer :: frazilSalinityTendency @@ -261,6 +262,7 @@ subroutine ocn_frazil_forcing_active_tracers(meshPool, tracersPool, forcingPool, call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(forcingPool, 'frazilTemperatureTendency', frazilTemperatureTendency) call mpas_pool_get_array(forcingPool, 'frazilSalinityTendency', frazilSalinityTendency) @@ -272,7 +274,7 @@ subroutine ocn_frazil_forcing_active_tracers(meshPool, tracersPool, forcingPool, !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) activeTracersTend(indexTemperature,k,iCell) = activeTracersTend(indexTemperature,k,iCell) & + frazilTemperatureTendency(k,iCell) activeTracersTend(indexSalinity,k,iCell) = activeTracersTend(indexSalinity,k,iCell) + frazilSalinityTendency(k,iCell) @@ -367,7 +369,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP real (kind=RKIND), pointer, dimension(:,:) :: layerThickness real (kind=RKIND), pointer, dimension(:,:,:) :: activeTracers - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell integer, pointer :: indexTemperature !< index in tracers array for temperature integer, pointer :: indexSalinity !< index in tracers array for salinity integer :: kBottomFrazil ! k index where testing for frazil begins @@ -393,6 +395,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) ! get mesh information + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) ! get arrays @@ -451,7 +454,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP ! find deepest level where frazil can be created kBottomFrazil=maxLevelCell(iCell) - do k=maxLevelCell(iCell), 1, -1 + do k=maxLevelCell(iCell), minLevelCell(iCell), -1 ! add the ssh so frazil can form below land ice, where the ssh is depressed if (-zMid(k,iCell) < -ssh(iCell) + config_frazil_maximum_depth) then kBottomFrazil=k @@ -461,7 +464,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP ! find minimum temperature between 1:kBottomFrazil columnTemperatureMin = 1.0e30_RKIND - do k=1,kBottomFrazil + do k=minLevelCell(iCell),kBottomFrazil if ( activeTracers(indexTemperature,k,iCell) < columnTemperatureMin ) then columnTemperatureMin = activeTracers(indexTemperature,k,iCell) end if @@ -475,7 +478,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP sumNewThicknessWeightedSaltContent = 0.0_RKIND ! loop from maximum depth of frazil creation to surface - do k = kBottomFrazil, 1, -1 + do k = kBottomFrazil, minLevelCell(iCell), -1 ! get freezing temperature oceanFreezingTemperature = ocn_freezing_temperature(salinity=activeTracers(indexSalinity,k,iCell), & @@ -570,7 +573,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP endif ! if (freezingEnergy < 0) - enddo ! do k=kBottom,1-1 + enddo ! do k=kBottom,minLevelCell,-1 ! accumulate frazil mass to column total ! note: the accumulatedFrazilIceMass (at both time levels) is reset to zero after being sent to the coupler From 57f09f9fe9949dd0b813c9f552890412be307e7e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 10:26:51 -0600 Subject: [PATCH 12/46] minLevel* changes to ocn_tracer_surface_flux_to_tend --- .../shared/mpas_ocn_tracer_surface_flux_to_tend.F | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_surface_flux_to_tend.F b/src/core_ocean/shared/mpas_ocn_tracer_surface_flux_to_tend.F index 4960fe8405..29bb3e46b9 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_surface_flux_to_tend.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_surface_flux_to_tend.F @@ -123,7 +123,7 @@ subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbso integer :: iCell, k, iTracer, nTracers, nCells integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND) :: remainingFlux @@ -137,6 +137,7 @@ subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbso call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) nTracers = size(tend, dim=1) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) nCells = nCellsArray( 1 ) @@ -145,7 +146,7 @@ subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbso !$omp do schedule(runtime) private(remainingFlux, k, iTracer) do iCell = 1, nCells remainingFlux = 1.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) remainingFlux = remainingFlux - fractionAbsorbed(k, iCell) do iTracer = 1, nTracers @@ -174,7 +175,7 @@ subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbso !$omp do schedule(runtime) private(remainingFlux, k, iTracer) do iCell = 1, nCells remainingFlux = 1.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) remainingFlux = remainingFlux - fractionAbsorbedRunoff(k, iCell) do iTracer = 1, nTracers From 82a35c48ac117e30ebcfef337beaa7d0095a809b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 15:45:59 -0600 Subject: [PATCH 13/46] minLevel* changes to ocn_tracer_nonlocalflux --- src/core_ocean/shared/mpas_ocn_tracer_nonlocalflux.F | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_nonlocalflux.F b/src/core_ocean/shared/mpas_ocn_tracer_nonlocalflux.F index 3d16139d40..668650ef07 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_nonlocalflux.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_nonlocalflux.F @@ -112,7 +112,7 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace integer :: iCell, k, iTracer, nTracers, nCells integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND) :: fluxTopOfCell, fluxBottomOfCell err = 0 @@ -125,6 +125,7 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) nTracers = size(tend, dim=1) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) nCells = nCellsArray( 1 ) @@ -132,7 +133,7 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace !$omp parallel !$omp do schedule(runtime) private(k, iTracer, fluxTopOfCell, fluxBottomOfCell) do iCell = 1, nCells - do k = 2, maxLevelCell(iCell)-1 + do k = minLevelCell(iCell)+1, maxLevelCell(iCell)-1 ! NOTE: at the moment, all tracers are based on the flux-profile used for temperature, i.e. vertNonLocalFlux(1,:,:) do iTracer = 1, nTracers @@ -151,7 +152,7 @@ subroutine ocn_tracer_nonlocalflux_tend(meshPool, vertNonLocalFlux, surfaceTrace end do ! enforce boundary conditions at top of column - k = 1 + k = minLevelCell(iCell) do iTracer = 1, nTracers fluxTopOfCell = 0.0_RKIND fluxBottomOfCell = surfaceTracerFlux(iTracer, iCell) * vertNonLocalFlux(1, k+1, iCell) From 44d8292de6c4f58c5049282a619b749cf138a8ee Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 16:34:14 -0600 Subject: [PATCH 14/46] minLevel* changes to ocn_diags,ocn_compute_kpp_rhs --- src/core_ocean/shared/mpas_ocn_diagnostics.F | 12 ++++++------ src/core_ocean/shared/mpas_ocn_vmix.F | 14 ++++++++------ 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics.F b/src/core_ocean/shared/mpas_ocn_diagnostics.F index 06101a3312..66ad9b923b 100644 --- a/src/core_ocean/shared/mpas_ocn_diagnostics.F +++ b/src/core_ocean/shared/mpas_ocn_diagnostics.F @@ -2016,8 +2016,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, timeLe ! ! Compute fraction of thickness flux that is in the top model layer - fracAbsorbed = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient, -100.0_RKIND) ) - fracAbsorbedRunoff = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient_runoff, & + fracAbsorbed = 1.0_RKIND - exp( max(-layerThickness(minLevelCell(iCell), iCell) / config_flux_attenuation_coefficient, -100.0_RKIND) ) + fracAbsorbedRunoff = 1.0_RKIND - exp( max(-layerThickness(minLevelCell(iCell), iCell) / config_flux_attenuation_coefficient_runoff, & -100.0_RKIND) ) ! Store the total tracer flux below in nonLocalSurfaceTemperatureFlux for use in the CVMix nonlocal ! transport code. This includes tracer forcing due to thickness @@ -2029,12 +2029,12 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, timeLe - fracAbsorbedRunoff * activeTracersSurfaceFluxRunoff(indexTempFlux, iCell) nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) = activeTracersSurfaceFlux(indexSaltFlux,iCell) & - - fracAbsorbed * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,1,iCell) & - - fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * activeTracers(indexSaltFlux,1,iCell) + - fracAbsorbed * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,minLevelCell(iCell),iCell) & + - fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * activeTracers(indexSaltFlux,minLevelCell(iCell),iCell) - surfaceBuoyancyForcing(iCell) = thermalExpansionCoeff (1,iCell) & + surfaceBuoyancyForcing(iCell) = thermalExpansionCoeff (minLevelCell(iCell),iCell) & * nonLocalSurfaceTracerFlux(indexTempFlux,iCell) & - - salineContractionCoeff(1,iCell) * nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) + - salineContractionCoeff(minLevelCell(iCell),iCell) * nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) ! at this point, surfaceBuoyancyForcing has units of m/s ! change into units of m^2/s^3 (which can be thought of as the flux of buoyancy, units of buoyancy * velocity ) diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index a859541f6a..6c326b6bd9 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -975,9 +975,9 @@ subroutine ocn_tracer_vmix_tend_implicit(meshPool, dt, vertDiffTopOfCell, layerT enddo if ( config_cvmix_kpp_nonlocal_with_implicit_mix ) then - call ocn_compute_kpp_rhs(tracers(:,:,iCell), rhs(:,:), dt, N, num_tracers, & + call ocn_compute_kpp_rhs(tracers(:,:,iCell), rhs(:,:), dt, Nsurf, N, num_tracers, & layerThickness(:,iCell), vertNonLocalFlux(:,:,iCell), & - tracerGroupSurfaceFlux(:,iCell)) !mlc-arguments may need to be changed + tracerGroupSurfaceFlux(:,iCell)) else rhs(:,:) = tracers(:,:,iCell) endif @@ -1376,24 +1376,26 @@ end subroutine tridiagonal_solve_mult!}}} ! !----------------------------------------------------------------------- -subroutine ocn_compute_kpp_rhs(tracers, rhs, dt, maxLevelCell, nTracers, & +subroutine ocn_compute_kpp_rhs(tracers, rhs, dt, minLevelCell, maxLevelCell, nTracers, & layerThickness, vertNonLocalFlux, tracerGroupSurfaceFlux)!{{{ real (kind=RKIND), intent(in) :: dt real (kind=RKIND), dimension(:,:), intent(in) :: vertNonLocalFlux, tracers real (kind=RKIND), dimension(:), intent(in) :: layerThickness, tracerGroupSurfaceFlux real (kind=RKIND), dimension(:,:), intent(out) :: rhs - integer, intent(in) :: maxLevelCell, nTracers + integer, intent(in) :: minLevelCell, maxLevelCell, nTracers integer :: iTracer, k - do k=2,maxLevelCell-1 + rhs(:,1:minLevelCell-1) = 0.0_RKIND + + do k=minLevelCell+1,maxLevelCell-1 do iTracer=1,nTracers rhs(iTracer, k) = tracers(iTracer,k) + dt * tracerGroupSurfaceFlux(iTracer) * & (vertNonLocalFlux(1,k) - vertNonLocalFlux(1,k+1)) / layerThickness(k) enddo enddo - k=1 + k=minLevelCell do iTracer=1,nTracers rhs(iTracer, k) = tracers(iTracer,k) + dt * tracerGroupSurfaceFlux(iTracer) * & (-vertNonLocalFlux(1,k+1) )/ layerThickness(k) From aa47979eb0896be1f7f77d7b753a7802c6945207 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 24 Mar 2021 09:06:51 -0600 Subject: [PATCH 15/46] Assign unused values in ocn_vmix --- src/core_ocean/shared/mpas_ocn_vmix.F | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index 6c326b6bd9..5e906db6dd 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -956,13 +956,14 @@ subroutine ocn_tracer_vmix_tend_implicit(meshPool, dt, vertDiffTopOfCell, layerT N = maxLevelCell(iCell) ! A is lower diagonal term - A(Nsurf)=0.0_RKIND + A(1:Nsurf)=0.0_RKIND do k = Nsurf+1, N A(k) = -2.0_RKIND*dt*vertDiffTopOfCell(k,iCell) & / (layerThickness(k-1,iCell) + layerThickness(k,iCell)) / layerThickness(k,iCell) enddo ! C is upper diagonal term + C(1:Nsurf-1)=0.0_RKIND do k = Nsurf, N-1 C(k) = -2.0_RKIND*dt*vertDiffTopOfCell(k+1,iCell) & / (layerThickness(k,iCell) + layerThickness(k+1,iCell)) / layerThickness(k,iCell) @@ -970,6 +971,7 @@ subroutine ocn_tracer_vmix_tend_implicit(meshPool, dt, vertDiffTopOfCell, layerT C(N) = 0.0_RKIND ! B is diagonal term + B(1:Nsurf-1)=0.0_RKIND do k = Nsurf, N B(k) = 1.0_RKIND - A(k) - C(k) enddo From e423bd6ee3576e58c27a15401e27ff3ec48728ed Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:32:26 -0600 Subject: [PATCH 16/46] minLevel* changes to ocn_tracer_advection_std --- .../shared/mpas_ocn_tracer_advection.F | 4 ++-- .../shared/mpas_ocn_tracer_advection_std.F | 20 +++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection.F b/src/core_ocean/shared/mpas_ocn_tracer_advection.F index ec49b4d7f2..899e5876a5 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection.F @@ -152,8 +152,8 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & else call ocn_tracer_advection_std_tend(tracers, advCoefs, advCoefs3rd, & nAdvCellsForEdge, advCellsForEdge, normalThicknessFlux, w, layerThickness, & - layerThickness, dt, meshPool, tend, maxLevelCell, maxLevelEdgeTop, & - highOrderAdvectionMask, edgeSignOnCell) + layerThickness, dt, meshPool, tend, minLevelCell, maxLevelCell, & + minLevelEdgeBot, maxLevelEdgeTop, highOrderAdvectionMask, edgeSignOnCell) endif call mpas_timer_stop("tracer adv") diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F index a34a88e1f4..c6aee55e4c 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F @@ -65,7 +65,7 @@ module ocn_tracer_advection_std !----------------------------------------------------------------------- subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdvCellsForEdge, advCellsForEdge, &!{{{ normalThicknessFlux, w, layerThickness, verticalCellSize, dt, meshPool, & - tend, maxLevelCell, maxLevelEdgeTop, & + tend, minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop, & highOrderAdvectionMask, edgeSignOnCell) real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: current tracer values @@ -80,8 +80,8 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv real (kind=RKIND), intent(in) :: dt !< Input: Timestep type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency - integer, dimension(:), pointer :: maxLevelCell !< Input: Index to max level at cell center - integer, dimension(:), pointer :: maxLevelEdgeTop !< Input: Index to max level at edge with non-land cells on both sides + integer, dimension(:), pointer :: minLevelCell, maxLevelCell !< Input: Index to min,max level at cell center + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop !< Input: Index to min,max level at edge with non-land cells on both sides integer, dimension(:,:), pointer :: highOrderAdvectionMask !< Input: Mask for high order advection integer, dimension(:, :), pointer :: edgeSignOnCell !< Input: Sign for flux from edge on each cell. @@ -145,12 +145,12 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Compute the high order vertical flux. Also determine bounds on tracer_cur. !$omp do schedule(runtime) private(k, verticalWeightK, verticalWeightKm1) do iCell = 1, nCells - k = max(1, min(maxLevelCell(iCell), 2)) + k = max(minLevelCell(iCell), min(maxLevelCell(iCell), minLevelCell(iCell)+1)) verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell)) - do k=3,maxLevelCell(iCell)-1 + do k=minLevelCell(iCell)+2,maxLevelCell(iCell)-1 select case (vertOrder) case (vertOrder4) high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), & @@ -166,7 +166,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv end select ! vertOrder end do - k = max(1, maxLevelCell(iCell)) + k = max(minLevelCell(iCell), maxLevelCell(iCell)) verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell)) @@ -180,7 +180,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv cell2 = cellsOnEdge(2, iEdge) ! Compute 2nd order fluxes where needed. - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5_RKIND) & * normalThicknessFlux(k, iEdge) @@ -191,7 +191,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Compute 3rd or 4th fluxes where requested. do i = 1, nAdvCellsForEdge(iEdge) iCell = advCellsForEdge(i,iEdge) - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef3rdOrder & * sign(1.0_RKIND,normalThicknessFlux(k,iEdge))*adv_coefs_3rd(i,iEdge)) @@ -208,7 +208,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv invAreaCell1 = 1.0_RKIND / areaCell(iCell) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) & * invAreaCell1 end do @@ -219,7 +219,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Accumulate the scaled high order vertical tendencies. !$omp do schedule(runtime) private(k) do iCell = 1, nCellsSolve - do k = 1,maxLevelCell(iCell) + do k = minLevelCell(iCell),maxLevelCell(iCell) tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) & - high_order_vert_flux(k, iCell)) end do ! k loop From 5112d947690dfab8a2316a1efd2729a13ffb33f5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 16:25:24 -0600 Subject: [PATCH 17/46] minLevel* changes to ocn_wetting_drying --- .../shared/mpas_ocn_wetting_drying.F | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_wetting_drying.F b/src/core_ocean/shared/mpas_ocn_wetting_drying.F index 7344ecfced..7c4436ed5f 100644 --- a/src/core_ocean/shared/mpas_ocn_wetting_drying.F +++ b/src/core_ocean/shared/mpas_ocn_wetting_drying.F @@ -108,7 +108,7 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{ !----------------------------------------------------------------- type (mpas_pool_type), pointer :: statePool, meshPool, tendPool - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND), dimension(:), pointer :: sshSubcycleNew real (kind=RKIND), dimension(:), pointer :: bottomDepth integer, pointer :: nCellsSolve @@ -134,6 +134,7 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{ call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, timeLevel=1) call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, timeLevel=2) call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) call mpas_pool_get_array(statePool, 'sshSubcycle', sshSubcycleNew, 2) @@ -148,7 +149,7 @@ subroutine ocn_wetting_drying_verify( block , minHeight, err)!{{{ ! check to make sure that there is no layer that is too dry minThickness = +1.0E34 do iCell = 1, nCellsSolve - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) ! use ssh as a proxy too for baroclinic mode if (trim(config_time_integrator) == 'split_explicit' .or. trim(config_time_integrator) == 'semi_implicit') then layerThick = min(layerThicknessNew(k, iCell), (sshSubcycleNew(iCell)+bottomDepth(iCell))/maxLevelCell(iCell)) @@ -243,6 +244,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + integer, dimension(:), pointer :: minLevelEdgeTop integer, dimension(:), pointer :: maxLevelEdgeTop integer, dimension(:), pointer :: maxLevelEdgeBot integer, pointer :: nEdges @@ -261,6 +263,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying call mpas_pool_get_array(provisStatePool, 'layerThickness', layerThicknessProvis, 1) call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) + call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', minLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) @@ -281,7 +284,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeBot(iEdge) + do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) < eps) then ! prevent spurious flux for close to zero values normalTransportVelocity(k, iEdge) = 0.0_RKIND @@ -370,8 +373,8 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, laye integer :: iEdge, iCell, k, i integer, pointer :: nVertLevels, nCells integer, dimension(:), pointer :: nEdgesOnCell - integer, dimension(:), pointer :: maxLevelCell - integer, dimension(:), pointer :: maxLevelEdgeBot + integer, dimension(:), pointer :: minLevelCell, maxLevelCell + integer, dimension(:), pointer :: minLevelEdgeTop, maxLevelEdgeBot integer, dimension(:,:), pointer :: edgesOnCell integer, dimension(:,:), pointer :: edgeSignOnCell @@ -387,7 +390,9 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, laye call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', minLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -399,12 +404,12 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, laye do iCell = 1, nCells invAreaCell = 1.0_RKIND / areaCell(iCell) ! can switch with maxLevelEdgeBot(iEdge) - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) divOutFlux = 0.0_RKIND layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell)) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge)) then + if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then ! only consider divergence flux leaving the cell if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then divOutFlux = divOutFlux + normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) & @@ -421,7 +426,7 @@ subroutine ocn_wetting_drying_wettingVelocity(meshPool, layerThicknessEdge, laye do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - if (k <= maxLevelEdgeBot(iEdge)) then + if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then ! each outgoing velocity is penalized (but not the incoming, wetting velocities) ! square the fractional term to make values near zero go to zero much quicker (to prevent threshold from being hit) From 17fbfeb2c2f7da1a2ac722e167831f7fef65a18f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 16:06:47 -0600 Subject: [PATCH 18/46] minLevel* changes to ocn_tidal_forcing --- src/core_ocean/shared/mpas_ocn_tidal_forcing.F | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F index 54cb42267e..e73fbae12a 100644 --- a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F @@ -107,7 +107,7 @@ subroutine ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, layerThickne integer :: iCell, k, nCells integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real (kind=RKIND), dimension(:,:), pointer :: tidalLayerThicknessTendency err = 0 @@ -117,6 +117,7 @@ subroutine ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, layerThickne call mpas_timer_start("tidal thickness tendency") call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(forcingPool, 'tidalLayerThicknessTendency', & tidalLayerThicknessTendency) @@ -128,7 +129,7 @@ subroutine ocn_tidal_forcing_layer_thickness(meshPool, forcingPool, layerThickne !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessTend(k,iCell) = layerThicknessTend(k,iCell) + & tidalLayerThicknessTendency(k,iCell) @@ -200,7 +201,7 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo integer :: iCell, k, nCells integer, dimension(:), pointer :: nCellsArray - integer, pointer, dimension(:) :: maxLevelCell + integer, pointer, dimension(:) :: minLevelCell, maxLevelCell integer, pointer :: nVertLevels real (kind=RKIND) :: dt @@ -225,6 +226,7 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue) call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) @@ -256,16 +258,16 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo ! compute total depth for relative thickness contribution totalDepth = 0.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) totalDepth = totalDepth + layerThickness(k,iCell) end do - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tidalLayerThicknessTendency(:,iCell) = 0.0_RKIND end do if (trim(config_tidal_forcing_type) == 'thickness_source') then ! distribute tidal forcing tendency fractionally over water column - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tidalLayerThicknessTendency(k,iCell) = tidalInputMask(iCell) / config_use_tidal_forcing_tau & * (layerThickness(k,iCell)/totalDepth) * (tidalHeight - ssh(iCell)) end do From 08e2af127990288459d3ac8a03418426511d9ec2 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:34:21 -0600 Subject: [PATCH 19/46] minLevel* changes to ocn_tracer_hmix_del2 --- src/core_ocean/shared/mpas_ocn_tracer_hmix_del2.F | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_hmix_del2.F b/src/core_ocean/shared/mpas_ocn_tracer_hmix_del2.F index 1cb65b179e..e8ac99ae3d 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_hmix_del2.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_hmix_del2.F @@ -119,7 +119,7 @@ subroutine ocn_tracer_hmix_del2_tend(meshPool, layerThicknessEdge, tracers, tend integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray - integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, nEdgesOnCell integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell real (kind=RKIND) :: invAreaCell @@ -137,6 +137,7 @@ subroutine ocn_tracer_hmix_del2_tend(meshPool, layerThicknessEdge, tracers, tend call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) num_tracers = size(tracers, dim=1) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'areaCell', areaCell) @@ -166,7 +167,7 @@ subroutine ocn_tracer_hmix_del2_tend(meshPool, layerThicknessEdge, tracers, tend r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) do iTracer = 1, num_tracers ! \kappa_2 \nabla \phi on edge tracer_turb_flux = tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1) From b2013a212651005dfa7330d1558459a7573eaef9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 13:35:56 -0600 Subject: [PATCH 20/46] minLevel* changes to ocn_tracer_hmix_del4 --- src/core_ocean/shared/mpas_ocn_tracer_hmix_del4.F | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_hmix_del4.F b/src/core_ocean/shared/mpas_ocn_tracer_hmix_del4.F index 8e22887ee5..6ece3a72f6 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_hmix_del4.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_hmix_del4.F @@ -120,7 +120,7 @@ subroutine ocn_tracer_hmix_del4_tend(meshPool, layerThicknessEdge, tracers, tend integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray, nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, minLevelCell, maxLevelCell, nEdgesOnCell integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge, edgesOnCell, edgeSignOnCell real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2 @@ -149,6 +149,8 @@ subroutine ocn_tracer_hmix_del4_tend(meshPool, layerThicknessEdge, tracers, tend call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) num_tracers = size(tracers, dim=1) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) @@ -181,7 +183,7 @@ subroutine ocn_tracer_hmix_del4_tend(meshPool, layerThicknessEdge, tracers, tend cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) do iTracer = 1, num_tracers * edgeMask(k, iEdge) r_tmp1 = invdcEdge * layerThicknessEdge(k, iEdge) * tracers(iTracer, k, cell1) @@ -211,7 +213,7 @@ subroutine ocn_tracer_hmix_del4_tend(meshPool, layerThicknessEdge, tracers, tend invdcEdge = meshScalingDel4(iEdge) * dvEdge(iEdge) * eddyDiff4 / dcEdge(iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) do iTracer = 1, num_tracers * edgeMask(k, iEdge) tracer_turb_flux = (delsq_tracer(iTracer, k, cell2) - delsq_tracer(iTracer, k, cell1)) From 5c4f3168eea56226b8d3f9f9940b05818944c210 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 14:53:49 -0600 Subject: [PATCH 21/46] minLevel* changes to ocn_tracer_hmix_redi --- .../shared/mpas_ocn_tracer_hmix_redi.F | 32 +++++++++++++------ 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F index 05f4649f90..e93b92086f 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F @@ -133,7 +133,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray, nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell, maxLevelCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, nEdgesOnCell, minLevelCell, maxLevelCell integer, dimension(:, :), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell, cellsOnCell real(kind=RKIND) :: tempTracer, invAreaCell1, invAreaCell2, invAreaCell, areaEdge @@ -151,6 +151,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) nTracers = size(tracers, dim=1) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) @@ -221,7 +223,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg r_tmp = dvEdge(iEdge)/dcEdge(iEdge) coef = dvEdge(iEdge) - k = 1 + k = minLevelEdgeBot(iEdge) kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* & 0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + & @@ -239,12 +241,22 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg /(zMid(k, cell1) - zMid(k + 1, cell1)) & + slopeTriadDown(k, 2, iEdge)*(tracers(iTr, k, cell2) - tracers(iTr, k + 1, cell2)) & /(zMid(k, cell2) - zMid(k + 1, cell2))) + if (minLevelCell(cell1) < minLevelEdgeBot(iEdge)) then + flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & + layerThicknessEdge(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 1, iEdge)* & + (tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1))/(zMid(k - 1, cell1) - zMid(k, cell1)) + endif + if (minLevelCell(cell2) < minLevelEdgeBot(iEdge)) then + flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & + layerThicknessEdge(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 2, iEdge)* & + (tracers(iTr, k - 1, cell2) - tracers(iTr, k, cell2))/(zMid(k - 1, cell2) - zMid(k, cell2)) + endif redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2*invAreaCell redi_term2_edge(iTr, k, iEdge) = -flux_term2 end do - do k = 2, maxLevelEdgeTop(iEdge) - 1 + do k = minLevelEdgeBot(iEdge) + 1, maxLevelEdgeTop(iEdge) - 1 kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))*0.25_RKIND* & @@ -324,15 +336,15 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg end do ! nEdgesOnCell(iCell) ! impose no-flux boundary conditions at top and bottom of column - fluxRediZTop(:, 1) = 0.0_RKIND + fluxRediZTop(:, 1:minLevelCell(iCell)) = 0.0_RKIND fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND - redi_term3_topOfCell(:, 1, iCell) = 0.0_RKIND - do k = 2, maxLevelCell(iCell) + redi_term3_topOfCell(:, 1:minLevelCell(iCell), iCell) = 0.0_RKIND + do k = minLevelCell(iCell) + 1, maxLevelCell(iCell) redi_term3_topOfCell(:, k, iCell) = fluxRediZTop(:, k) * invAreaCell end do redi_term3_topOfCell(:, maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, nTracers ! Add tendency for Term 3: d/dz ( h S grad phi) = ( S grad phi) fluxes ! 2.0 in next line is because a dot product on a C-grid @@ -367,7 +379,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg !$omp parallel !$omp do schedule(runtime) private(k, iTr, tempTracer, iEdge) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tempTracer = tracers(iTr, k, iCell) + dt*(redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + & redi_term3(iTr, k, iCell))/layerThickness(k, iCell) @@ -396,7 +408,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg !$omp do schedule(runtime) private(invAreaCell, k, iTr, i, iEdge) do iCell = 1, nCells invAreaCell = 1.0_RKIND/areaCell(iCell) - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tend(iTr, k, iCell) = tend(iTr, k, iCell) + redi_term1(iTr, k, iCell) tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* & @@ -409,7 +421,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg do i = 1, nEdgesOnCell(iCell) if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle iEdge = edgesOnCell(i, iCell) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) do iTr = 1, ntracers tend(iTr, k, iCell) = tend(iTr, k, iCell) + edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge) end do From 7eda6ec8e7566b36385d2d06e1770d73b355073f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 26 Mar 2021 17:44:08 -0600 Subject: [PATCH 22/46] fix to Redi needed to ensure bfb results in ecosys tracers --- src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F index e93b92086f..eeb8a12faf 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F @@ -224,6 +224,11 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg coef = dvEdge(iEdge) k = minLevelEdgeBot(iEdge) + kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & + RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* & + 0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + & + RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2)) + k = 1 kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))* & 0.25_RKIND*(RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + & From 842ba3524fe08932b5495a8a7cf410bae31cd927 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 15:40:04 -0600 Subject: [PATCH 23/46] minLevel* changes to ocn_vel_hmix_del4 --- src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F index 94ed413f50..9eba8dc2de 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F @@ -150,7 +150,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !Compute Laplacian of velocity #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc present(cellsOnEdge, verticesOnEdge, minLevelEdgeBot, maxLevelEdgeTop, & !$acc dvEdge, dcEdge, del2u, div, relVort) & !$acc private(k, cell1, cell2, vertex1, vertex2, & !$acc dcEdgeInv, dvEdgeInv) @@ -171,7 +171,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) dvEdgeInv = 1.0_RKIND / max(dvEdge(iEdge), 0.25_RKIND*dcEdge(iEdge)) - do k=1,maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge) ! Compute \nabla^2 u = \nabla divergence + k \times \nabla relativeVorticity del2u(k,iEdge) = (div(k,cell2) - div(k,cell1))*dcEdgeInv & -(relVort(k,vertex2) - relVort(k,vertex1))*& @@ -188,7 +188,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ ! Compute del2 of relative vorticity #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(edgesOnVertex, maxLevelVertexTop, dcEdge, & + !$acc present(edgesOnVertex, minLevelVertexBot, maxLevelVertexTop, dcEdge, & !$acc areaTriangle, edgeSignOnVertex, del2u, & !$acc del2relVort) & !$acc private(i, k, iEdge, areaTriInv) @@ -201,7 +201,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ areaTriInv = 1.0_RKIND/areaTriangle(iVertex) do i = 1, vertexDegree iEdge = edgesOnVertex(i, iVertex) - do k = 1, maxLevelVertexTop(iVertex) + do k = minLevelVertexBot(iVertex), maxLevelVertexTop(iVertex) del2relVort(k,iVertex) = del2relVort(k,iVertex) & + edgeSignOnVertex(i,iVertex) & *dcEdge(iEdge)*del2u(k,iEdge)*areaTriInv @@ -218,7 +218,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ ! Compute del2 of divergence #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(nEdgesOnCell, edgesOnCell, maxLevelCell, dvEdge,& + !$acc present(nEdgesOnCell, edgesOnCell, minLevelCell, maxLevelCell, dvEdge,& !$acc edgeSignOnCell, areaCell, del2u, del2div) & !$acc private(i, k, iEdge, areaCellInv) #else @@ -230,7 +230,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ areaCellInv = 1.0_RKIND / areaCell(iCell) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) del2div(k,iCell) = del2div(k,iCell) - & edgeSignOnCell(i,iCell)*dvEdge(iEdge) & *del2u(k,iEdge)*areaCellInv @@ -248,7 +248,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc present(cellsOnEdge, verticesOnEdge, minLevelEdgeBot, maxLevelEdgeTop, & !$acc dcEdge, dvEdge, meshScalingDel4, edgeMask, & !$acc del2div, del2relVort, tend) & !$acc private(k, cell1, cell2, vertex1, vertex2, & @@ -269,7 +269,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ dvEdgeInv = 1.0_RKIND / dvEdge(iEdge) visc4 = viscDel4*meshScalingDel4(iEdge) - do k=1, maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) uDiff = divFactor*(del2div(k,cell2) - del2div(k,cell1))* & dcEdgeInv & - (del2relVort(k,vertex2) - del2relVort(k,vertex1))* & From ad58e5b22e1e2019e77a1d7650d0eebc064ee6ec Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 15:44:02 -0600 Subject: [PATCH 24/46] minLevel* changes to ocn_vel_hmix_leith --- src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F index 467ab15f19..ed005162b0 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F @@ -139,7 +139,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc present(cellsOnEdge, verticesOnEdge, minLevelEdgeBot, maxLevelEdgeTop, & !$acc dcEdge, dvEdge, meshScalingDel2, div, relVort, & !$acc tend, edgeMask) & !$acc private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & @@ -161,7 +161,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ visc2tmp = (leithParam*dxLeith*meshScalingDel2(iEdge)/pii)**3 - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) ! Here -( relativeVorticity(k,vertex2) - ! relativeVorticity(k,vertex1) ) / dvEdge(iEdge) From 1ad68816e46b704a5f200bc373b06a6ddec3f584 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 15:53:59 -0600 Subject: [PATCH 25/46] minLevel* changes to ocn_vel_vmix_tend_rayleigh,... ocn_vel_vmix_tend_spatially_variable, ocn_vel_vmix_tend_spatially_variable_manning --- src/core_ocean/shared/mpas_ocn_vmix.F | 88 ++++++++++++++++----------- 1 file changed, 53 insertions(+), 35 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index 5e906db6dd..3bfd4e01c8 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -243,11 +243,11 @@ subroutine ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, N, nEdges + integer :: iEdge, k, cell1, cell2, Nsurf, N, nEdges integer, pointer :: nVertLevels integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop integer, dimension(:,:), pointer :: cellsOnEdge @@ -261,18 +261,19 @@ subroutine ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) nEdges = nEdgesArray( 1 ) allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) - A(1)=0.0_RKIND !$omp parallel !$omp do schedule(runtime) & - !$omp private(N, cell1, cell2, edgeThicknessTotal, k, A, B, C, velTemp) + !$omp private(Nsurf, N, cell1, cell2, edgeThicknessTotal, k, A, B, C, velTemp) do iEdge = 1, nEdges + Nsurf = minLevelEdgeBot(iEdge) N = maxLevelEdgeTop(iEdge) if (N .gt. 0) then @@ -282,30 +283,33 @@ subroutine ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) edgeThicknessTotal = 0.0_RKIND - do k = 1, N + do k = Nsurf, N layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) edgeThicknessTotal = edgeThicknessTotal + layerThicknessEdge(k,iEdge) end do ! A is lower diagonal term - do k = 2, N + A(1:Nsurf)=0.0_RKIND + do k = Nsurf+1, N A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! C is upper diagonal term - do k = 1, N-1 + C(1:Nsurf-1)=0.0_RKIND + do k = Nsurf, N-1 C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! B is diagonal term - B(1) = 1.0_RKIND - C(1) & + B(1:Nsurf-1)=0.0_RKIND + B(Nsurf) = 1.0_RKIND - C(Nsurf) & ! added Rayleigh terms + dt*rayleighDampingCoef/((1.0_RKIND - rayleighDepthVariable) + rayleighDepthVariable*edgeThicknessTotal) - do k = 2, N-1 + do k = Nsurf, N-1 B(k) = 1.0_RKIND - A(k) - C(k) & ! added Rayleigh terms + dt*(rayleighDampingCoef/((1.0_RKIND - rayleighDepthVariable) + rayleighDepthVariable*edgeThicknessTotal)) @@ -319,9 +323,10 @@ subroutine ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, + dt*(rayleighBottomDampingCoef + & rayleighDampingCoef /((1.0_RKIND - rayleighDepthVariable) + rayleighDepthVariable*edgeThicknessTotal)) - call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + call tridiagonal_solve(A(Nsurf+1:N),B(Nsurf:N),C(Nsurf:N-1),normalVelocity(Nsurf:N,iEdge),velTemp(Nsurf:N),N-Nsurf+1) - normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(1:Nsurf-1,iEdge) = 0.0_RKIND + normalVelocity(Nsurf:N,iEdge) = velTemp(Nsurf:N) normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND end if @@ -424,7 +429,7 @@ subroutine ocn_vel_vmix_tend_implicit(meshPool, dt, kineticEnergyCell, vertViscT allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) !$omp parallel - !$omp do schedule(runtime) private(N, cell1, cell2, A, B, C, velTemp, k) + !$omp do schedule(runtime) private(Nsurf, N, cell1, cell2, A, B, C, velTemp, k) do iEdge = 1, nEdges N = maxLevelEdgeTop(iEdge) Nsurf = minLevelEdgeBot(iEdge) @@ -554,12 +559,12 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, N, nEdges + integer :: iEdge, k, cell1, cell2, N, Nsurf, nEdges integer, pointer :: nVertLevels real (kind=RKIND) :: implicitCd integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop integer, dimension(:,:), pointer :: cellsOnEdge @@ -572,16 +577,18 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) nEdges = nEdgesArray( 1 ) allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) - A(1)=0.0_RKIND - !$omp do schedule(runtime) + !$omp parallel + !$omp do schedule(runtime) private(Nsurf, N, cell1, cell2, A, B, C, velTemp, k) do iEdge = 1, nEdges + Nsurf = minLevelEdgeBot(iEdge) N = maxLevelEdgeTop(iEdge) if (N .gt. 0) then @@ -590,7 +597,7 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d ! so recompute layerThicknessEdge here. cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k = 1, N + do k = Nsurf, N layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) end do @@ -598,22 +605,25 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d implicitCd = 0.5_RKIND*(bottomDrag(cell1) + bottomDrag(cell2)) ! A is lower diagonal term - do k = 2, N + A(1:Nsurf)=0.0_RKIND + do k = Nsurf+1, N A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! C is upper diagonal term - do k = 1, N-1 + C(1:Nsurf-1)=0.0_RKIND + do k = Nsurf, N-1 C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! B is diagonal term - B(1) = 1.0_RKIND - C(1) - do k = 2, N-1 + B(1:Nsurf-1)=0.0_RKIND + B(Nsurf) = 1.0_RKIND - C(Nsurf) + do k = Nsurf+1, N-1 B(k) = 1.0_RKIND - A(k) - C(k) enddo @@ -623,9 +633,10 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d B(N) = 1.0_RKIND - A(N) + dt*implicitCd & * sqrt(kineticEnergyCell(N,cell1) + kineticEnergyCell(N,cell2)) / layerThicknessEdge(N,iEdge) - call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + call tridiagonal_solve(A(Nsurf+1:N),B(Nsurf:N),C(Nsurf:N-1),normalVelocity(Nsurf:N,iEdge),velTemp(Nsurf:N),N-Nsurf+1) - normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(1:Nsurf-1,iEdge) = 0.0_RKIND + normalVelocity(Nsurf:N,iEdge) = velTemp(Nsurf:N) normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND end if @@ -717,12 +728,12 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, N, nEdges + integer :: iEdge, k, cell1, cell2, N, Nsurf, nEdges integer, pointer :: nVertLevels real (kind=RKIND) :: implicitCd integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop integer, dimension(:,:), pointer :: cellsOnEdge @@ -748,6 +759,7 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) @@ -761,7 +773,6 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc nCells = nCellsSolve allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) - A(1)=0.0_RKIND ! Compute bottomDrag (Manning roughness) induced by vegetation if (config_use_vegetation_drag .AND. config_use_vegetation_manning_equation) then @@ -794,8 +805,10 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc enddo endif - !$omp do schedule(runtime) + !$omp parallel + !$omp do schedule(runtime) private(Nsurf, N, A, B, C, rhs, velTemp, k) do iEdge = 1, nEdges + Nsurf = minLevelEdgeBot(iEdge) N = maxLevelEdgeTop(iEdge) if (N .gt. 0) then @@ -804,7 +817,7 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc ! so recompute layerThicknessEdge here. cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k = 1, N + do k = Nsurf, N layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) end do @@ -818,22 +831,25 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc endif ! A is lower diagonal term - do k = 2, N + A(1:Nsurf)=0.0_RKIND + do k = Nsurf+1, N A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! C is upper diagonal term - do k = 1, N-1 + C(1:Nsurf-1)=0.0_RKIND + do k = Nsurf, N-1 C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & / layerThicknessEdge(k,iEdge) enddo ! B is diagonal term - B(1) = 1.0_RKIND - C(1) - do k = 2, N-1 + B(1:Nsurf-1)=0.0_RKIND + B(Nsurf) = 1.0_RKIND - C(Nsurf) + do k = Nsurf+1, N-1 B(k) = 1.0_RKIND - A(k) - C(k) enddo @@ -843,14 +859,16 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc B(N) = 1.0_RKIND - A(N) + dt*implicitCd & * sqrt(kineticEnergyCell(N,cell1) + kineticEnergyCell(N,cell2)) / layerThicknessEdge(N,iEdge) - call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + call tridiagonal_solve(A(Nsurf+1:N),B(Nsurf:N),C(Nsurf:N-1),normalVelocity(Nsurf:N,iEdge),velTemp(Nsurf:N),N-Nsurf+1) - normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(1:Nsurf-1,iEdge) = 0.0_RKIND + normalVelocity(Nsurf:N,iEdge) = velTemp(Nsurf:N) normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND end if end do !$omp end do + !$omp end parallel deallocate(A,B,C,velTemp) @@ -949,7 +967,7 @@ subroutine ocn_tracer_vmix_tend_implicit(meshPool, dt, vertDiffTopOfCell, layerT call mpas_timer_start('vmix tracers tend imp loop', .false.) !$omp parallel - !$omp do schedule(runtime) private(N, A, B, C, rhs, tracersTemp, k) + !$omp do schedule(runtime) private(Nsurf, N, A, B, C, rhs, tracersTemp, k) do iCell = 1, nCells ! Compute A(k), B(k), C(k) for tracers Nsurf = minLevelCell(iCell) From 398d22f48dae4b1ff6d7ea7cb2f85533f05fd97a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Mar 2021 15:37:14 -0600 Subject: [PATCH 26/46] minLevel* changes to ocn_gm --- src/core_ocean/shared/mpas_ocn_gm.F | 73 +++++++++++++++-------------- 1 file changed, 39 insertions(+), 34 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_gm.F b/src/core_ocean/shared/mpas_ocn_gm.F index 7f6266302e..8ed920a0e9 100644 --- a/src/core_ocean/shared/mpas_ocn_gm.F +++ b/src/core_ocean/shared/mpas_ocn_gm.F @@ -124,7 +124,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND), dimension(:), pointer :: & ssh, fEdge real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, tridiagA, tridiagB, tridiagC, rightHandSide - integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, minLevelCell, maxLevelCell, nEdgesOnCell integer, dimension(:, :), pointer :: cellsOnEdge, edgesOnCell integer :: i, k, iEdge, cell1, cell2, iCell, N, iter, iCellSelf, maxLocation real(kind=RKIND) :: h1, h2, areaEdge, c, BruntVaisalaFreqTopEdge, rtmp @@ -167,6 +167,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) @@ -236,7 +238,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule(runtime) private(maxLocation, k, BruntVaisalaFreqTopEdge, maxN) do iCell = 1, nCells - k = min(maxLevelCell(iCell) - 1, max(1, indMLD(iCell))) + k = min(maxLevelCell(iCell) - 1, max(minLevelCell(iCell), indMLD(iCell))) maxN = max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND) do while (BruntVaisalaFreqTop(k + 1, iCell) > maxN .and. k < maxLevelCell(iCell) - 1) k = k + 1 @@ -259,10 +261,10 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule (runtime) private(zMLD, k) do iCell = 1, nCells - k = max(1, indMLD(iCell)) + k = max(minLevelCell(iCell), indMLD(iCell)) zMLD = ssh(iCell) - zMid(k, iCell) - do k = 1, indMLD(iCell) + do k = minLevelCell(iCell), indMLD(iCell) RediKappaSfcTaper(k, iCell) = abs((ssh(iCell) - zMid(k, iCell))/(zMLD)) end do end do @@ -281,7 +283,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & k33(1:maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND k33Norm(1:maxLevelCell(iCell) + 1) = epsGM ! prep dz, dTdz and dSdz for this column - do k = 2, maxLevelCell(iCell) + do k = minLevelCell(iCell)+1, maxLevelCell(iCell) dzTop(k) = 0.5_RKIND*(layerThickness(k - 1, iCell) + layerThickness(k, iCell)) dTdzTop(k) = (activeTracers(indexTemperature, k - 1, iCell) & - activeTracers(indexTemperature, k, iCell)) & @@ -290,9 +292,9 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & - activeTracers(indexSalinity, k, iCell)) & /dzTop(k) end do - dzTop(1) = -1e-15_RKIND - dTdzTop(1) = -1e-15_RKIND - dSdzTop(1) = -1e-15_RKIND + dzTop(1:minLevelCell(iCell)) = -1e-15_RKIND + dTdzTop(1:minLevelCell(iCell)) = -1e-15_RKIND + dSdzTop(1:minLevelCell(iCell)) = -1e-15_RKIND dzTop(maxLevelCell(iCell) + 1) = -1e-15_RKIND dTdzTop(maxLevelCell(iCell) + 1) = -1e-15_RKIND dSdzTop(maxLevelCell(iCell) + 1) = -1e-15_RKIND @@ -309,7 +311,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & dcEdgeInv = 1.0_RKIND/dcEdge(iEdge) areaEdge = dcEdge(iEdge)*dvEdge(iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) drhoDT = -thermExpCoeff(k, iCell) drhoDS = salineContractCoeff(k, iCell) dTdx = (activeTracers(indexTemperature, k, cell2) & @@ -365,7 +367,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & slopeTaperDown*sfcTaperDown*slopeTriadDown(k, iCellSelf, iEdge) ! Griffies 1998 eqn 34 - if (k > 1) then + if (k > minLevelCell(iCell)) then k33(k, iCell) = k33(k, iCell) + & areaEdge*dzTop(k)*slopeTriadUp(k, iCellSelf, iEdge)**2 k33Norm(k) = k33Norm(k) + areaEdge*dzTop(k) @@ -381,10 +383,10 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & end do ! nEdgesOnCell(iCell) ! Normalize k33 - do k = 2, maxLevelCell(iCell) + do k = minLevelCell(iCell)+1, maxLevelCell(iCell) k33(k, iCell) = k33(k, iCell)/k33Norm(k)*RediKappaScaling(k, iCell) end do - k33(1, iCell) = 0.0_RKIND + k33(1:minLevelCell(iCell), iCell) = 0.0_RKIND k33(maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND end do ! iCell !$omp end do @@ -420,7 +422,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule(runtime) private(k, rtmp) do iCell = 1, nCells - do k = 2, maxLevelCell(iCell) + do k = minLevelCell(iCell)+1, maxLevelCell(iCell) rtmp = (displacedDensity(k-1,iCell) - density(k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell)) dDensityDzTopOfCell(k,iCell) = min(rtmp, -epsGM) end do @@ -428,7 +430,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & ! Approximation of dDensityDzTopOfCell on the top and bottom interfaces through the idea of having ! ghost cells above the top and below the bottom layers of the same depths and density. ! Essentially, this enforces the boundary condition (d density)/dz = 0 at the top and bottom. - dDensityDzTopOfCell(1,iCell) = 0.0_RKIND + dDensityDzTopOfCell(1:minLevelCell(iCell),iCell) = 0.0_RKIND dDensityDzTopOfCell(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND end do !$omp end do @@ -440,7 +442,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule(runtime) private(k, cell1, cell2) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeTop(iEdge)+1 + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1 cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) dDensityDzTopOfEdge(k,iEdge) = 0.5_RKIND * (dDensityDzTopOfCell(k,cell1) + dDensityDzTopOfCell(k,cell2)) @@ -458,7 +460,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k=1,maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge) gradDensityEdge(k,iEdge) = (density(k,cell2) - density(k,cell1)) / dcEdge(iEdge) gradZMidEdge(k,iEdge) = (zMid(k,cell2) - zMid(k,cell1)) / dcEdge(iEdge) end do @@ -469,7 +471,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & do iEdge = 1, nEdges ! The interpolation can only be carried out on non-boundary edges if (maxLevelEdgeTop(iEdge) .GE. 1) then - do k = 2, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) h1 = layerThickEdge(k-1,iEdge) h2 = layerThickEdge(k,iEdge) ! Using second-order interpolation below @@ -481,9 +483,9 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & ! Approximation of values on the top and bottom interfaces through the idea of having ghost cells ! above the top and below the bottom layers of the same depths and density. - gradDensityTopOfEdge(1,iEdge) = gradDensityEdge(1,iEdge) + gradDensityTopOfEdge(minLevelEdgeBot(iEdge),iEdge) = gradDensityEdge(minLevelEdgeBot(iEdge),iEdge) gradDensityTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradDensityEdge(maxLevelEdgeTop(iEdge),iEdge) - gradZMidTopOfEdge(1,iEdge) = gradZMidEdge(1,iEdge) + gradZMidTopOfEdge(minLevelEdgeBot(iEdge),iEdge) = gradZMidEdge(minLevelEdgeBot(iEdge),iEdge) gradZMidTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradZMidEdge(maxLevelEdgeTop(iEdge),iEdge) end if end do @@ -492,7 +494,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges if (maxLevelEdgeTop(iEdge) .GE. 1) then - do k = 1, maxLevelEdgeTop(iEdge)+1 + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1 gradDensityConstZTopOfEdge(k,iEdge) = gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & * gradZMidTopOfEdge(k,iEdge) end do @@ -508,19 +510,19 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule(runtime) private(k, BruntVaisalaFreqTopEdge, maxN, cell1, cell2) do iEdge=1,nEdges - k=1 + k=minLevelEdgeBot(iEdge) cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k,cell1) + & BruntVaisalaFreqTop(k,cell2)) maxN = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) - do k = 2, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k,cell1) + & BruntVaisalaFreqTop(k,cell2)) maxN = max(maxN,max(BruntVaisalaFreqTopEdge, 0.0_RKIND)) enddo - do k = 2, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k,cell1) + & BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) @@ -546,7 +548,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & cell2 = cellsOnEdge(2, iEdge) sumN2 = 0.0 - do k = 2, maxLevelEdgeTop(iEdge)-1 + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)-1 lt1 = 0.5_RKIND*(layerThickness(k,cell1) + layerThickness(k-1,cell1)) lt2 = 0.5_RKIND*(layerThickness(k,cell2) + layerThickness(k-1,cell2)) @@ -580,7 +582,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp do schedule(runtime) private(k, cell1, cell2, sumN2, ltSum, sumRi, countN2, & !$omp zEdge, RiTopOfEdge, BruntVaisalaFreqTopEdge, lt1, lt2) do iEdge = 1, nEdges - k=1 + k=minLevelEdgeBot(iEdge) sumN2 = 0.0_RKIND ltSum = 0.0_RKIND sumRi = 0.0_RKIND @@ -588,7 +590,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) - k = 2 + k = minLevelEdgeBot(iEdge)+1 zEdge = -layerThickEdge(k-1,iEdge) do while(zEdge > -config_GM_Visbeck_max_depth .and. k < maxLevelEdgeTop(iEdge)) lt1 = 0.5_RKIND*(layerThickness(k,cell1) + layerThickness(k-1,cell1)) @@ -637,7 +639,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & Lr = min(cGMphaseSpeed(iEdge) / abs(fEdge(iEdge)), & sqrt(cGMphaseSpeed(iEdge) / (2.0_RKIND*betaEdge(iEdge)))) - do k=2,maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge)+1,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND*(max(BruntVaisalaFreqTop(k,cell1),0.0_RKIND) + & max(BruntVaisalaFreqTop(k,cell2),0.0_RKIND)) shearEdgeInv = (0.5_RKIND*(layerThickEdge(k-1,iEdge) + layerThickEdge(k,iEdge)))**2.0 & @@ -673,7 +675,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & ! Construct the tridiagonal matrix if (maxLevelEdgeTop(iEdge) .GE. 3) then ! First row - k = 2 + k = minLevelEdgeBot(iEdge)+1 BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k, cell1) + BruntVaisalaFreqTop(k, cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) @@ -685,7 +687,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & *gradDensityConstZTopOfEdge(k,iEdge) ! Second to next to the last rows - do k = 3, maxLevelEdgeTop(iEdge) - 1 + do k = minLevelEdgeBot(iEdge) + 2, maxLevelEdgeTop(iEdge) - 1 BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k, cell1) + BruntVaisalaFreqTop(k, cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagA(k - 2) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdge(k - 1, iEdge) & @@ -709,11 +711,14 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw & *gradDensityConstZTopOfEdge(k,iEdge) ! Total number of rows - N = maxLevelEdgeTop(iEdge) - 1 + N = maxLevelEdgeTop(iEdge) - minLevelEdgeBot(iEdge) ! Call the tridiagonal solver - call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, & - gmStreamFuncTopOfEdge(2:maxLevelEdgeTop(iEdge), iEdge), N) + call tridiagonal_solve(tridiagA(minLevelEdgeBot(iEdge):maxLevelEdgeTop(iEdge)), & + tridiagB(minLevelEdgeBot(iEdge):maxLevelEdgeTop(iEdge)), & + tridiagC(minLevelEdgeBot(iEdge):maxLevelEdgeTop(iEdge)), & + rightHandSide(minLevelEdgeBot(iEdge):maxLevelEdgeTop(iEdge)), & + gmStreamFuncTopOfEdge(minLevelEdgeBot(iEdge)+1:maxLevelEdgeTop(iEdge), iEdge), N) end if end do !!$omp end do @@ -723,7 +728,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) normalGMBolusVelocity(k, iEdge) = gmResolutionTaper(iEdge) * (gmStreamFuncTopOfEdge(k, iEdge) & - gmStreamFuncTopOfEdge(k + 1, iEdge)) /layerThickEdge(k, iEdge) end do @@ -743,7 +748,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & areaEdge = 0.25_RKIND*dcEdge(iEdge)*dvEdge(iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) rtmp = 0.5_RKIND*(gmStreamFuncTopOfEdge(k, iEdge) + gmStreamFuncTopOfEdge(k + 1, iEdge))*areaEdge gmStreamFuncTopOfCell(k, iCell) = gmStreamFuncTopOfCell(k, iCell) + gmResolutionTaper(iEdge) * rtmp end do From 8d7ae8c9e202823cbd74ca7ea75c220bf1036e76 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Mar 2021 14:38:59 -0600 Subject: [PATCH 27/46] Minlevelcell changes to ocn_init_isomip_plus --- src/core_ocean/mode_init/mpas_ocn_init_isomip_plus.F | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/mode_init/mpas_ocn_init_isomip_plus.F b/src/core_ocean/mode_init/mpas_ocn_init_isomip_plus.F index 58a59b3e45..0a07c213bd 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_isomip_plus.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_isomip_plus.F @@ -104,7 +104,7 @@ subroutine ocn_init_setup_isomip_plus(domain, iErr)!{{{ integer, pointer :: nCells, nVertLevels integer, pointer :: index_temperature, index_salinity, index_tracer1 - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell real(kind=RKIND), dimension(:), pointer :: refBottomDepth, & fCell, fEdge, fVertex, effectiveDensityInLandIce, xCell, & @@ -186,6 +186,7 @@ subroutine ocn_init_setup_isomip_plus(domain, iErr)!{{{ call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) @@ -194,7 +195,7 @@ subroutine ocn_init_setup_isomip_plus(domain, iErr)!{{{ call mpas_pool_get_dimension(tracersPool, 'index_salinity', index_salinity) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) frac = (0.0_RKIND-zMid(k, iCell))/(0.0_RKIND-config_isomip_plus_max_bottom_depth) activeTracers(index_temperature, k, iCell) = (1.0_RKIND - frac)*config_isomip_plus_init_top_temp & + frac*config_isomip_plus_init_bot_temp From a7a03377569c7b7e27dbe4a196f6f0fc5e47f11a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Mar 2021 09:12:18 -0600 Subject: [PATCH 28/46] Fixup minLevel* changes to ocn_tendency --- src/core_ocean/shared/mpas_ocn_tendency.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 85818bd987..c0b64d7d6b 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -1037,7 +1037,7 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, meshPool, timeL call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', maxLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) From 51b0c4988a0e14e8119d8f3fc9754ede381d7da6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Mar 2021 14:05:04 -0600 Subject: [PATCH 29/46] minLevel* changes to ocn_time_integration_si --- .../mpas_ocn_time_integration_si.F | 56 ++++++++++--------- 1 file changed, 31 insertions(+), 25 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F index ee7acfb9a7..67fdf81ec4 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F @@ -157,7 +157,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ integer, dimension(:), pointer :: nCellsArray, nEdgesArray ! Mesh array pointers - integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, nEdgesOnEdge, nEdgesOnCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop, nEdgesOnEdge, nEdgesOnCell integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask, edgesOnEdge integer, dimension(:,:), pointer :: edgesOnCell, edgeSignOnCell @@ -492,6 +492,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_subpool(tendPool, 'tracersTend', tracersTendPool) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) @@ -518,7 +519,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ cell2 = cellsOnEdge(2,iEdge) uTemp = 0.0_RKIND ! could put this after with uTemp(maxleveledgetop+1:nvertlevels)=0 - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) ! normalBaroclinicVelocityNew = normalBaroclinicVelocityOld + dt*(-f*normalBaroclinicVelocityPerp ! + T(u*,w*,p*) + g*grad(SSH*) ) @@ -533,17 +534,17 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! thicknessSum is initialized outside the loop because on land boundaries ! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a ! nonzero value to avoid a NaN. - normalThicknessFluxSum = layerThickEdge(1,iEdge) * uTemp(1) - thicknessSum = layerThickEdge(1,iEdge) + normalThicknessFluxSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge) * uTemp(minLevelEdgeBot(iEdge)) + thicknessSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge) - do k = 2, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) normalThicknessFluxSum = normalThicknessFluxSum + layerThickEdge(k,iEdge) * uTemp(k) thicknessSum = thicknessSum + layerThickEdge(k,iEdge) enddo barotropicForcing(iEdge) = split * normalThicknessFluxSum / thicknessSum / dt - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) ! These two steps are together here: !{\bf u}'_{k,n+1} = {\bf u}'_{k,n} - \Delta t {\overline {\bf G}} !{\bf u}'_{k,n+1/2} = \frac{1}{2}\left({\bf u}^{'}_{k,n} +{\bf u}'_{k,n+1}\right) @@ -2436,6 +2437,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'normalBarotropicVelocity', normalBarotropicVelocityNew, 2) call mpas_pool_get_array(statePool, 'normalBaroclinicVelocity', normalBaroclinicVelocityNew, 2) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) @@ -2469,10 +2471,10 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! thicknessSum is initialized outside the loop because on land boundaries ! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a ! nonzero value to avoid a NaN. - normalThicknessFluxSum = layerThickEdge(1,iEdge) * uTemp(1) - thicknessSum = layerThickEdge(1,iEdge) + normalThicknessFluxSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge) * uTemp(minLevelEdgeBot(iEdge)) + thicknessSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge) - do k = 2, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) normalThicknessFluxSum = normalThicknessFluxSum + layerThickEdge(k,iEdge) * uTemp(k) thicknessSum = thicknessSum + layerThickEdge(k,iEdge) enddo @@ -2617,8 +2619,10 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool) call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(tracersPool, 'activeTracers', tracersGroupCur, 1) @@ -2664,7 +2668,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp do schedule(runtime) private(k, temp_h, temp, i) do iCell = 1, nCells ! sshNew is a pointer, defined above. - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) ! this is h_{n+1} temp_h = layerThicknessCur(k,iCell) + dt * layerThicknessTend(k,iCell) @@ -2689,7 +2693,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k, temp) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) ! h^{hf}_{n+1} was computed in Stage 1 @@ -2741,7 +2745,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) ! this is h_{n+1} layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell) + dt * layerThicknessTend(k,iCell) end do @@ -2753,7 +2757,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k= 1, maxLevelEdgeTop(iEdge) + do k= minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) activeTracerHorizontalAdvectionEdgeFlux(:,k,iEdge) = & activeTracerHorizontalAdvectionEdgeFlux(:,k,iEdge) / & layerThickEdge(k,iEdge) @@ -2763,7 +2767,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k= 1, maxLevelCell(iCell) + do k= minLevelCell(iCell), maxLevelCell(iCell) activeTracerHorizontalAdvectionTendency(:,k,iCell) = & activeTracerHorizontalAdvectionTendency(:,k,iCell) / & layerThicknessNew(k,iCell) @@ -2805,7 +2809,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersGroupNew(:,k,iCell) = (tracersGroupCur(:,k,iCell) * layerThicknessCur(k,iCell) + dt & * tracersGroupTend(:,k,iCell) ) / layerThicknessNew(k,iCell) end do @@ -2818,7 +2822,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersGroupNew(indexSalinity,k,iCell) = max(0.001_RKIND, tracersGroupNew(indexSalinity,k,iCell)) end do end do @@ -2834,7 +2838,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) ! h^{hf}_{n+1} was computed in Stage 1 @@ -2859,7 +2863,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) normalVelocityNew(k,iEdge) = normalBarotropicVelocityNew(iEdge) + 2 * normalBaroclinicVelocityNew(k,iEdge) & - normalBaroclinicVelocityCur(k,iEdge) end do @@ -3071,7 +3075,7 @@ subroutine ocn_time_integration_si_init(domain)!{{{ type (dm_info) :: dminfo integer :: iTracer, cell, cell1, cell2 - integer, dimension(:), pointer :: maxLevelEdgeTop,maxLevelCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, minLevelCell, maxLevelCell integer, dimension(:,:), pointer :: cellsOnEdge real (kind=RKIND) :: normalThicknessFluxSum, layerThicknessSum, layerThickEdge1 real (kind=RKIND), dimension(:), pointer :: refBottomDepth, normalBarotropicVelocity @@ -3140,7 +3144,9 @@ subroutine ocn_time_integration_si_init(domain)!{{{ call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'areaCell', areaCell) call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) @@ -3167,11 +3173,11 @@ subroutine ocn_time_integration_si_init(domain)!{{{ ! thicknessSum is initialized outside the loop because on land boundaries ! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a ! nonzero value to avoid a NaN. - layerThickEdge1 = 0.5_RKIND*( layerThickness(1,cell1) + layerThickness(1,cell2) ) - normalThicknessFluxSum = layerThickEdge1 * normalVelocity(1,iEdge) + layerThickEdge1 = 0.5_RKIND*( layerThickness(minLevelCell(cell1),cell1) + layerThickness(minLevelCell(cell2),cell2) ) + normalThicknessFluxSum = layerThickEdge1 * normalVelocity(minLevelEdgeBot(iEdge),iEdge) layerThicknessSum = layerThickEdge1 - do k=2, maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) ! ocn_diagnostic_solve has not yet been called, so compute hEdge ! just for this edge. layerThickEdge1 = 0.5_RKIND*( layerThickness(k,cell1) + layerThickness(k,cell2) ) @@ -3184,7 +3190,7 @@ subroutine ocn_time_integration_si_init(domain)!{{{ normalBarotropicVelocity(iEdge) = normalThicknessFluxSum / layerThicknessSum ! normalBaroclinicVelocity(k,iEdge) = normalVelocity(k,iEdge) - normalBarotropicVelocity(iEdge) - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) normalBaroclinicVelocity(k,iEdge) = normalVelocity(k,iEdge) - normalBarotropicVelocity(iEdge) enddo @@ -3272,12 +3278,12 @@ subroutine ocn_time_integration_si_init(domain)!{{{ k = maxLevelCell(iCell) zTop(k:nVertLevels,iCell) = -bottomDepth(iCell) + layerThickness(k,iCell) - do k = maxLevelCell(iCell)-1, 1, -1 + do k = maxLevelCell(iCell)-1, minLevelCell(iCell), -1 zTop(k,iCell) = zTop(k+1,iCell) + layerThickness(k ,iCell) end do ! copy zTop(1,iCell) into sea-surface height array - sshCur(iCell) = zTop(1,iCell) + sshCur(iCell) = zTop(minLevelCell(iCell),iCell) end do tmp1 = minval(sshCur) From 34b45b933ac79fcfc24175a24721f30320224c60 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Mar 2021 09:12:18 -0600 Subject: [PATCH 30/46] fixup ocn_tend_freq_filtered_thickness --- src/core_ocean/shared/mpas_ocn_tendency.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index c0b64d7d6b..6c6e580b9e 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -1037,7 +1037,6 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, meshPool, timeL call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) @@ -1072,7 +1071,8 @@ subroutine ocn_tend_freq_filtered_thickness(tendPool, statePool, meshPool, timeL do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + ! Note: changing loop limit to minLevelEdgeBot introduces non-bit-for-bit changes + do k = 1, maxLevelEdgeTop(iEdge) flux = layerThickEdge(k, iEdge) * normalVelocity(k, iEdge) * dvEdge(iEdge) * edgeSignOnCell(i, iCell) * invAreaCell div_hu(k) = div_hu(k) - flux div_hu_btr = div_hu_btr - flux From ca9acf2fc23dcbdb8ea2168408da4e90d9ef4ac3 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 29 Mar 2021 18:12:37 -0600 Subject: [PATCH 31/46] Remaining minLevel* changes to cvmix --- src/core_ocean/shared/mpas_ocn_vmix_cvmix.F | 88 +++++++++++---------- 1 file changed, 45 insertions(+), 43 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F index 1591c88786..cb04e70443 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F @@ -364,7 +364,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JC084iC05p02503 ! The standard atan formulation is approximated following equation (13) of Girones et al (2012) ! https://doi.org/10.1109/MSP.2012.2219677 - do k = 1,maxLevelCell(iCell) + do k = minLevelCell(iCell),maxLevelCell(iCell) x = (abs(cvmix_variables % zw_iface(k)) - config_cvmix_BryanLewis_transitionDepth) / & config_cvmix_BryanLewis_transitionWidth vertDiffTopOfCell(k, iCell) = vertDiffTopOfCell(k, iCell) + config_cvmix_BryanLewis_bl1 + & @@ -375,11 +375,13 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim endif ! fill Ri - RiSmoothed(1:maxLevelCell(iCell)) = RiTopOfCell(1:maxLevelCell(iCell),iCell) + RiSmoothed(1:minLevelCell(iCell)-1) = RiTopOfCell(minLevelCell(iCell),iCell) + RiSmoothed(minLevelCell(iCell):maxLevelCell(iCell)) = RiTopOfCell(minLevelCell(iCell):maxLevelCell(iCell),iCell) RiSmoothed(maxLevelCell(iCell)+1) = RiSmoothed(maxLevelCell(iCell)) RiTemp(1:maxLevelCell(iCell)+1) = RiSmoothed(1:maxLevelCell(iCell)+1) ! Use a 1-2-1 filter to remove 2 dz noise in RiTopOfCell + ! Note: loop limits here must not be minLevelCell or non-bfb changes are introduced do nsmooth=1,config_cvmix_num_ri_smooth_loops do k=2,maxLevelCell(iCell) RiSmoothed(k) = (RiTemp(k-1) + 2.0_RKIND*RiTemp(k) + RiTemp(k+1)) / 4.0_RKIND @@ -387,7 +389,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim RiTemp(1:maxLevelCell(iCell)) = RiSmoothed(1:maxLevelCell(iCell)) enddo - cvmix_variables%ShearRichardson_iface => RiSmoothed(1:maxLevelCell(iCell)+1) + cvmix_variables%ShearRichardson_iface => RiSmoothed(minLevelCell(iCell):maxLevelCell(iCell)+1) endif @@ -396,17 +398,17 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim BruntVaisalaFreqTop(minLevelCell(iCell):maxLevelCell(iCell),iCell)) BVFSmoothed(1:minLevelCell(iCell)-1) = max(0.0_RKIND,BVFSmoothed(minLevelCell(iCell))) BVFSmoothed(maxLevelCell(iCell)+1) = max(0.0_RKIND,BVFSmoothed(maxLevelCell(iCell))) - cvmix_variables%SqrBuoyancyFreq_iface => BVFSmoothed(1:maxLevelCell(iCell)+1) + cvmix_variables%SqrBuoyancyFreq_iface => BVFSmoothed(minLevelCell(iCell):maxLevelCell(iCell)+1) ! fill the intent(in) KPP cvmix_variables % SurfaceFriction = surfaceFrictionVelocity(iCell) cvmix_variables % SurfaceBuoyancyForcing = surfaceBuoyancyForcing(iCell) - cvmix_variables % BulkRichardson_cntr => bulkRichardsonNumber(1:maxLevelCell(iCell), iCell) + cvmix_variables % BulkRichardson_cntr => bulkRichardsonNumber(minLevelCell(iCell):maxLevelCell(iCell), iCell) if (kpp_stage == 2) then if (config_use_cvmix_shear) then - do k = 1, maxLevelCell(iCell) + 1 + do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 cvmix_variables % Mdiff_iface(k) = 0.0_RKIND cvmix_variables % Tdiff_iface(k) = 0.0_RKIND end do @@ -418,13 +420,13 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! accounted for seperately. so remove bac kground from shear mixing values if(cvmixShearMixingChoice==cvmixShearMixingTypePP) then - do k = 1, maxLevelCell(iCell) + 1 + do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 cvmix_variables % Mdiff_iface(k) = cvmix_variables % Mdiff_iface(k) - config_cvmix_background_viscosity cvmix_variables % Tdiff_iface(k) = cvmix_variables % Tdiff_iface(k) - config_cvmix_background_diffusion end do endif - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) vertViscTopOfCell(k, iCell) = vertViscTopOfCell(k, iCell) + cvmix_variables % Mdiff_iface(k) vertDiffTopOfCell(k, iCell) = vertDiffTopOfCell(k, iCell) + cvmix_variables % Tdiff_iface(k) end do @@ -438,8 +440,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if (config_use_cvmix_fixed_boundary_layer) then cvmix_variables % BoundaryLayerDepth = config_cvmix_kpp_boundary_layer_depth cvmix_variables % kOBL_depth = cvmix_kpp_compute_kOBL_depth( & - zw_iface = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1), & - zt_cntr = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & + zw_iface = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + zt_cntr = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & OBL_depth = cvmix_variables % BoundaryLayerDepth ) else @@ -448,7 +450,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim do k=minLevelCell(iCell),maxLevelCell(iCell) Nsqr_iface(k) = BVFSmoothed(k) enddo - Nsqr_iface(1:minLevelCell(iCell)-1) = Nsqr_iface(minLevelCell(iCell)) + Nsqr_iface(minLevelCell(iCell):minLevelCell(iCell)-1) = Nsqr_iface(minLevelCell(iCell)) k=min(maxLevelCell(iCell)+1,nVertLevels) Nsqr_iface(k:maxLevelCell(iCell)+1) = Nsqr_iface(k-1) @@ -507,17 +509,17 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! compute the turbulent scales in order to compute the bulk Richardson number call cvmix_kpp_compute_turbulent_scales( & sigma_coord = config_cvmix_kpp_surface_layer_extent, & - OBL_depth = OBLDepths(1:maxLevelCell(iCell)), & - surf_buoy_force = interfaceForcings(1:maxLevelCell(iCell)), & + OBL_depth = OBLDepths(minLevelCell(iCell):maxLevelCell(iCell)), & + surf_buoy_force = interfaceForcings(minLevelCell(iCell):maxLevelCell(iCell)), & surf_fric_vel = cvmix_variables % SurfaceFriction, & - w_s = turbulentScalarVelocityScale(1:maxLevelCell(iCell))) + w_s = turbulentScalarVelocityScale(minLevelCell(iCell):maxLevelCell(iCell))) ! averaging over a surface layer assuming that BLdepth is cell bottom ! Build deltaVelocitySquared do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - deltaVelocitySquared(1:minLevelCell(iCell)) = 0.0_RKIND + deltaVelocitySquared(minLevelCell(iCell):minLevelCell(iCell)) = 0.0_RKIND normalVelocitySum(minLevelCell(iCell)) = normalVelocity(minLevelCell(iCell), iEdge)* & layerThickEdge(minLevelCell(iCell),iEdge) @@ -558,11 +560,11 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! call mpas_timer_stop('Bulk Richardson kIndexOBL loops') cvmix_variables % bulkRichardson_cntr(:) = cvmix_kpp_compute_bulk_Richardson( & - zt_cntr = cvmix_variables % zt_cntr(1:maxLevelCell(iCell)), & - delta_buoy_cntr = bulkRichardsonNumberBuoy(1:maxLevelCell(iCell),iCell), & - delta_Vsqr_cntr = bulkRichardsonNumberShear(1:maxLevelCell(iCell),iCell), & + zt_cntr = cvmix_variables % zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & + delta_buoy_cntr = bulkRichardsonNumberBuoy(minLevelCell(iCell):maxLevelCell(iCell),iCell), & + delta_Vsqr_cntr = bulkRichardsonNumberShear(minLevelCell(iCell):maxLevelCell(iCell),iCell), & ws_cntr = turbulentScalarVelocityScale(:), & - Nsqr_iface = Nsqr_iface(1:maxLevelCell(iCell)+1), & + Nsqr_iface = Nsqr_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & EFactor = langmuirEnhancementFactor, & LaSL = langmuirNumber, & bfsfc = cvmix_variables%SurfaceBuoyancyForcing, & @@ -572,11 +574,11 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! each level of bulk Richardson is computed as if OBL resided at bottom of that level call cvmix_kpp_compute_OBL_depth( & - Ri_bulk = bulkRichardsonNumber(1:maxLevelCell(iCell),iCell), & - zw_iface = cvmix_variables % zw_iface(1:maxLevelCell(iCell)+1), & + Ri_bulk = bulkRichardsonNumber(minLevelCell(iCell):maxLevelCell(iCell),iCell), & + zw_iface = cvmix_variables % zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & OBL_depth = cvmix_variables % BoundaryLayerDepth, & kOBL_depth = cvmix_variables % kOBL_depth, & - zt_cntr = cvmix_variables % zt_cntr(1:maxLevelCell(iCell)), & + zt_cntr = cvmix_variables % zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & surf_fric = cvmix_variables % SurfaceFriction, & surf_buoy = cvmix_variables % SurfaceBuoyancyForcing, & Coriolis = cvmix_variables % Coriolis) @@ -587,8 +589,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if(cvmix_variables % BoundaryLayerDepth .lt. layerThickness(minLevelCell(iCell),iCell)/2.0_RKIND) then cvmix_variables % BoundaryLayerDepth = layerThickness(minLevelCell(iCell),iCell)/2.0_RKIND cvmix_variables % kOBL_depth = cvmix_kpp_compute_kOBL_depth( & - zw_iface = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1),& - zt_cntr = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & + zw_iface = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1),& + zt_cntr = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & OBL_depth = cvmix_variables % BoundaryLayerDepth ) endif @@ -597,8 +599,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if(cvmix_variables % BoundaryLayerDepth .lt. configure_cvmix_kpp_minimum_OBL_under_sea_ice) then cvmix_variables % BoundaryLayerDepth = configure_cvmix_kpp_minimum_OBL_under_sea_ice cvmix_variables % kOBL_depth = cvmix_kpp_compute_kOBL_depth( & - zw_iface = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1),& - zt_cntr = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & + zw_iface = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1),& + zt_cntr = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & OBL_depth = cvmix_variables % BoundaryLayerDepth ) endif endif @@ -607,8 +609,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if(cvmix_variables % BoundaryLayerDepth .gt. abs(cvmix_variables%zt_cntr(maxLevelCell(iCell)))) then cvmix_variables % BoundaryLayerDepth = abs(cvmix_variables%zt_cntr(maxLevelCell(iCell))) cvmix_variables % kOBL_depth = cvmix_kpp_compute_kOBL_depth( & - zw_iface = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1), & - zt_cntr = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & + zw_iface = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + zt_cntr = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & OBL_depth = cvmix_variables % BoundaryLayerDepth ) endif @@ -618,7 +620,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if (kpp_stage == 2) then ! copy data into cvmix_variables - do k = 1, maxLevelCell(iCell) + 1 + do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 cvmix_variables % Mdiff_iface(k) = vertViscTopOfCell(k, iCell) cvmix_variables % Tdiff_iface(k) = vertDiffTopOfCell(k, iCell) end do @@ -628,8 +630,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim boundaryLayerDepth(iCell) = min(blTemp, abs(cvmix_variables%zt_cntr(maxLevelCell(iCell)))) cvmix_variables % kOBL_depth = cvmix_kpp_compute_kOBL_depth( & - zw_iface = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1), & - zt_cntr = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & + zw_iface = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + zt_cntr = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & OBL_depth = boundaryLayerDepth(iCell) ) ! update Langmuir enhancement factor @@ -646,18 +648,18 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! call mpas_timer_start('cvmix coeffs kpp', .false.) call cvmix_coeffs_kpp( & - Mdiff_out = cvmix_variables % Mdiff_iface(1:maxLevelCell(iCell)+1), & - Tdiff_out = cvmix_variables % Tdiff_iface(1:maxLevelCell(iCell)+1), & - Sdiff_out = cvmix_variables % Sdiff_iface(1:maxLevelCell(iCell)+1), & - zw = cvmix_variables%zw_iface(1:maxLevelCell(iCell)+1), & - zt = cvmix_variables%zt_cntr(1:maxLevelCell(iCell)), & - old_Mdiff = cvmix_variables%Mdiff_iface(1:maxLevelCell(iCell)+1), & - old_Tdiff = cvmix_variables%Tdiff_iface(1:maxLevelCell(iCell)+1), & - old_Sdiff = cvmix_variables%Sdiff_iface(1:maxLevelCell(iCell)+1), & + Mdiff_out = cvmix_variables % Mdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + Tdiff_out = cvmix_variables % Tdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + Sdiff_out = cvmix_variables % Sdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + zw = cvmix_variables%zw_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + zt = cvmix_variables%zt_cntr(minLevelCell(iCell):maxLevelCell(iCell)), & + old_Mdiff = cvmix_variables%Mdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + old_Tdiff = cvmix_variables%Tdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + old_Sdiff = cvmix_variables%Sdiff_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & OBL_depth = boundaryLayerDepth(iCell), & kOBL_depth = cvmix_variables%kOBL_depth, & - Tnonlocal = cvmix_variables%kpp_Tnonlocal_iface(1:maxLevelCell(iCell)+1), & - Snonlocal = cvmix_variables%kpp_Snonlocal_iface(1:maxLevelCell(iCell)+1), & + Tnonlocal = cvmix_variables%kpp_Tnonlocal_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & + Snonlocal = cvmix_variables%kpp_Snonlocal_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & surf_fric = cvmix_variables%SurfaceFriction, & surf_buoy = cvmix_variables%SurfaceBuoyancyForcing, & nlev = maxLevelCell(iCell), & @@ -668,7 +670,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim ! intent out of BoundaryLayerDepth is boundary layer depth measured in meters and vertical index indexBoundaryLayerDepth(iCell) = cvmix_variables % kOBL_depth - do k = 1, maxLevelCell(iCell) + 1 + do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 vertViscTopOfCell(k, iCell) = cvmix_variables % Mdiff_iface(k) vertDiffTopOfCell(k, iCell) = cvmix_variables % Tdiff_iface(k) end do @@ -685,7 +687,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if ( kpp_stage == 2) then ! call convective mixing scheme if (config_use_cvmix_convection) then - do k = 1, maxLevelCell(iCell) + 1 + do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 cvmix_variables % Mdiff_iface(k) = 0.0_RKIND cvmix_variables % Tdiff_iface(k) = 0.0_RKIND end do From 6d3770f40d3eee246ed046c044fa192f965bd1fc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 31 Mar 2021 10:43:41 -0600 Subject: [PATCH 32/46] Fixup minLevelCell changes to vmix: Define vertViscTopOfEdge from minLevelEdgeBot --- src/core_ocean/shared/mpas_ocn_vmix.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index 3bfd4e01c8..c86cde817a 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -1051,7 +1051,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, real (kind=RKIND), dimension(:,:), pointer :: nonLocalSurfaceTracerFlux, tracerGroupSurfaceFlux real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroup real (kind=RKIND), dimension(:,:,:), allocatable :: nonLocalFluxTend - integer, dimension(:), pointer :: maxLevelCell, minLevelCell, maxLevelEdgeTop, minLevelEdgeTop + integer, dimension(:), pointer :: maxLevelCell, minLevelCell, maxLevelEdgeTop, minLevelEdgeBot integer, dimension(:,:), pointer :: cellsOnEdge type (mpas_pool_iterator_type) :: groupItr @@ -1078,7 +1078,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'minLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', minLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) @@ -1104,7 +1104,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool, vertViscTopOfEdge(:, iEdge) = 0.0_RKIND cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k=minLevelEdgeTop(iEdge),maxLevelEdgeTop(iEdge) + do k=minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge) vertViscTopOfEdge(k,iEdge) = 0.5_RKIND*(vertViscTopOfCell(k,cell2)+vertViscTopOfCell(k,cell1)) end do end do From 0ba1a67d59fde7e09e5a4d284b1ff616163a3092 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 31 Mar 2021 21:10:21 -0600 Subject: [PATCH 33/46] Fixup minLevel* changes to vmix, rayleigh Co-authored-by: Matthew Turner --- src/core_ocean/shared/mpas_ocn_vmix.F | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index c86cde817a..7810b5e0c5 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -307,9 +307,8 @@ subroutine ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, ! B is diagonal term B(1:Nsurf-1)=0.0_RKIND B(Nsurf) = 1.0_RKIND - C(Nsurf) & - ! added Rayleigh terms + dt*rayleighDampingCoef/((1.0_RKIND - rayleighDepthVariable) + rayleighDepthVariable*edgeThicknessTotal) - do k = Nsurf, N-1 + do k = Nsurf+1, N-1 B(k) = 1.0_RKIND - A(k) - C(k) & ! added Rayleigh terms + dt*(rayleighDampingCoef/((1.0_RKIND - rayleighDepthVariable) + rayleighDepthVariable*edgeThicknessTotal)) @@ -585,8 +584,7 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, d allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) - !$omp parallel - !$omp do schedule(runtime) private(Nsurf, N, cell1, cell2, A, B, C, velTemp, k) + !$omp do schedule(runtime) do iEdge = 1, nEdges Nsurf = minLevelEdgeBot(iEdge) N = maxLevelEdgeTop(iEdge) @@ -805,8 +803,7 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc enddo endif - !$omp parallel - !$omp do schedule(runtime) private(Nsurf, N, A, B, C, rhs, velTemp, k) + !$omp do schedule(runtime) do iEdge = 1, nEdges Nsurf = minLevelEdgeBot(iEdge) N = maxLevelEdgeTop(iEdge) @@ -868,7 +865,6 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc end if end do !$omp end do - !$omp end parallel deallocate(A,B,C,velTemp) From 4349195a81cc6307f0fb614a73254bbbb3862ad7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 6 Apr 2021 14:55:24 -0600 Subject: [PATCH 34/46] Resolve merge conflict with freshwater flux terms --- src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F index ac9b07287d..ba9b687cab 100644 --- a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F @@ -546,19 +546,19 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer do iCell = 1, nCells ! Accumulate fluxes that use the surface temperature - rainTemperatureFlux(iCell) = rainFlux(iCell) * tracerGroup(index_temperature_flux,1,iCell) / rho_sw - evapTemperatureFlux(iCell) = evaporationFlux(iCell) * tracerGroup(index_temperature_flux,1,iCell) / rho_sw + rainTemperatureFlux(iCell) = rainFlux(iCell) * tracerGroup(index_temperature_flux,minLevelCell(iCell),iCell) / rho_sw + evapTemperatureFlux(iCell) = evaporationFlux(iCell) * tracerGroup(index_temperature_flux,minLevelCell(iCell),iCell) / rho_sw ! Runoff can only have a minimum temperature of 0.0C, since it is fresh water. tracersSurfaceFluxRunoff(index_temperature_flux,iCell) = riverRunoffFlux(iCell) & - * max(tracerGroup(index_temperature_flux,1,iCell), 0.0_RKIND) / rho_sw + * max(tracerGroup(index_temperature_flux,minLevelCell(iCell),iCell), 0.0_RKIND) / rho_sw ! Accumulate fluxes that use the freezing point seaIceFreshWaterTemperatureFlux(iCell) = seaIceFreshWaterFlux(iCell) * & - ocn_freezing_temperature( tracerGroup(index_salinity_flux, 1, iCell), pressure=0.0_RKIND, & + ocn_freezing_temperature( tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell), pressure=0.0_RKIND, & inLandIceCavity=.false.) / rho_sw icebergFreshWaterTemperatureFlux(iCell) = icebergFreshWaterFlux(iCell) * & - ocn_freezing_temperature( tracerGroup(index_salinity_flux, 1, iCell), pressure=0.0_RKIND, & + ocn_freezing_temperature( tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell), pressure=0.0_RKIND, & inLandIceCavity=.false.) / rho_sw totalFreshWaterTemperatureFlux(iCell) = rainTemperatureFlux(iCell) + evapTemperatureFlux(iCell) + & From 21464cc575f4eda71e2c81d3d64da74c8e388aac Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 10:02:37 -0600 Subject: [PATCH 35/46] Fixup freq-filtered diagnostics --- src/core_ocean/shared/mpas_ocn_diagnostics.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics.F b/src/core_ocean/shared/mpas_ocn_diagnostics.F index 66ad9b923b..af8957d7cb 100644 --- a/src/core_ocean/shared/mpas_ocn_diagnostics.F +++ b/src/core_ocean/shared/mpas_ocn_diagnostics.F @@ -1859,13 +1859,13 @@ subroutine ocn_filter_btr_mode_tend_vel(tendPool, statePool, meshPool, timeLevel tend_normalVelocity(maxLevelEdgeBot(iEdge),iEdge) thicknessSum = layerThickEdge(maxLevelEdgeBot(iEdge),iEdge) - do k = maxLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) normalThicknessFluxSum = normalThicknessFluxSum + layerThickEdge(k,iEdge) * tend_normalVelocity(k,iEdge) thicknessSum = thicknessSum + layerThickEdge(k,iEdge) enddo vertSum = normalThicknessFluxSum / thicknessSum - do k = maxLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) tend_normalVelocity(k,iEdge) = tend_normalVelocity(k,iEdge) - vertSum enddo enddo ! iEdge From 4e2dbee723e774f036f23da949cc7f1cfcd0fb61 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 10:15:50 -0600 Subject: [PATCH 36/46] Fixup gm k-limit conditionals --- src/core_ocean/shared/mpas_ocn_gm.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_gm.F b/src/core_ocean/shared/mpas_ocn_gm.F index 8ed920a0e9..6bf3407e9f 100644 --- a/src/core_ocean/shared/mpas_ocn_gm.F +++ b/src/core_ocean/shared/mpas_ocn_gm.F @@ -470,7 +470,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp do schedule(runtime) private(k, h1, h2) do iEdge = 1, nEdges ! The interpolation can only be carried out on non-boundary edges - if (maxLevelEdgeTop(iEdge) .GE. 1) then + if (maxLevelEdgeTop(iEdge) .GE. minLevelEdgeBot(iEdge)) then do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) h1 = layerThickEdge(k-1,iEdge) h2 = layerThickEdge(k,iEdge) @@ -493,7 +493,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - if (maxLevelEdgeTop(iEdge) .GE. 1) then + if (maxLevelEdgeTop(iEdge) .GE. minLevelEdgeBot(iEdge)) then do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1 gradDensityConstZTopOfEdge(k,iEdge) = gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & * gradZMidTopOfEdge(k,iEdge) @@ -672,8 +672,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & gmStreamFuncTopOfEdge(:, iEdge) = 0.0_RKIND - ! Construct the tridiagonal matrix - if (maxLevelEdgeTop(iEdge) .GE. 3) then + ! Construct the tridiagonal matrix if there are at least 3 active cells + if ((maxLevelEdgeTop(iEdge) - minLevelEdgeBot(iEdge)) .GE. 2) then ! First row k = minLevelEdgeBot(iEdge)+1 From 430ab933810ce95ccecc414c16f6c89f9701dea8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 10:23:37 -0600 Subject: [PATCH 37/46] Fixup ocn_tidal_forcing --- src/core_ocean/shared/mpas_ocn_tidal_forcing.F | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F index e73fbae12a..c25909b039 100644 --- a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F @@ -262,9 +262,7 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePoo totalDepth = totalDepth + layerThickness(k,iCell) end do - do k = minLevelCell(iCell), maxLevelCell(iCell) - tidalLayerThicknessTendency(:,iCell) = 0.0_RKIND - end do + tidalLayerThicknessTendency(:,iCell) = 0.0_RKIND if (trim(config_tidal_forcing_type) == 'thickness_source') then ! distribute tidal forcing tendency fractionally over water column do k = minLevelCell(iCell), maxLevelCell(iCell) From ca5f52d2445504df7d2b2e50857d00c0d05b87b0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 10:25:25 -0600 Subject: [PATCH 38/46] Fixup ocn_tracer_DMS --- src/core_ocean/shared/mpas_ocn_tracer_DMS.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_DMS.F b/src/core_ocean/shared/mpas_ocn_tracer_DMS.F index 0076c867fb..7285f9ab33 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_DMS.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_DMS.F @@ -172,10 +172,10 @@ subroutine ocn_tracer_DMS_compute(activeTracers, DMSTracers, nTracersDMS, ecosys numColumns = 1 column = 1 - iLevelSurface = minLevelCell(iCell) !DWJ 08/05/2016: This loop needs OpenMP added to it. do iCell=1,nCellsSolve + iLevelSurface = minLevelCell(iCell) DMS_input%number_of_active_levels(column) = maxLevelCell(iCell) - minLevelCell(iCell) + 1 DMS_forcing%ShortWaveFlux_surface(column) = shortWaveHeatFlux(iCell) From 318fb1c78cc70edbfe82aad8829fedbe251b970d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 13:16:36 -0600 Subject: [PATCH 39/46] Fixup cvmix --- src/core_ocean/shared/mpas_ocn_vmix_cvmix.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F index cb04e70443..bccd588920 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F @@ -519,7 +519,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - deltaVelocitySquared(minLevelCell(iCell):minLevelCell(iCell)) = 0.0_RKIND + deltaVelocitySquared(1:minLevelCell(iCell)) = 0.0_RKIND normalVelocitySum(minLevelCell(iCell)) = normalVelocity(minLevelCell(iCell), iEdge)* & layerThickEdge(minLevelCell(iCell),iEdge) @@ -662,8 +662,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim Snonlocal = cvmix_variables%kpp_Snonlocal_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & surf_fric = cvmix_variables%SurfaceFriction, & surf_buoy = cvmix_variables%SurfaceBuoyancyForcing, & - nlev = maxLevelCell(iCell), & - max_nlev = nVertLevels, & + nlev = maxLevelCell(iCell) - minLevelCell(iCell) + 2, & + max_nlev = maxLevelCell(iCell) - minLevelCell(iCell) + 2, & Langmuir_EFactor = langmuirEnhancementFactor ) ! call mpas_timer_stop('cvmix coeffs kpp') @@ -687,7 +687,7 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim if ( kpp_stage == 2) then ! call convective mixing scheme if (config_use_cvmix_convection) then - do k = minLevelCell(iCell), maxLevelCell(iCell) + 1 + do k = 1, maxLevelCell(iCell) + 1 cvmix_variables % Mdiff_iface(k) = 0.0_RKIND cvmix_variables % Tdiff_iface(k) = 0.0_RKIND end do From bcdcde32ceb34e3267e66194a345b684b063686b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 29 Mar 2021 18:37:41 -0600 Subject: [PATCH 40/46] minLevel* changes to RK4 timestepping --- .../mpas_ocn_time_integration_rk4.F | 40 +++++++++++-------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F index 781c5babff..f4e2d6f676 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F @@ -141,7 +141,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ ! Diagnostics Indices ! Mesh array pointers - integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop real (kind=RKIND), dimension(:), pointer :: bottomDepth ! Provis Array Pointers @@ -234,13 +234,15 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergenceCur, 1) call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergenceNew, 2) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) normalVelocityNew(k, iEdge) = normalVelocityCur(k, iEdge) end do end do @@ -248,7 +250,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessNew(k, iCell) = layerThicknessCur(k, iCell) end do end do @@ -267,7 +269,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells ! couple tracers to thickness - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersNew(:, k, iCell) = tracersCur(:, k, iCell) * layerThicknessCur(k, iCell) end do end do @@ -658,12 +660,12 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ if (tidalInputMask(iCell) == 1.0_RKIND) then ! compute total depth for relative thickness contribution totalDepth = 0.0_RKIND - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) totalDepth = totalDepth + restingThickness(k,iCell) end do ! only modify layer thicknesses on tidal boundary - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessNew(k, iCell) = tidalInputMask(iCell)*(tidalBCValue(iCell) + bottomDepth(iCell))*(restingThickness(k,iCell)/totalDepth) !(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell) ! generalized tappered assumption code end do @@ -719,8 +721,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells - surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(1, iCell) - surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(1, iCell) + surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(minLevelEdgeBot(iEdge), iCell) + surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(minLevelEdgeBot(iEdge), iCell) SSHGradient(indexSSHGradientZonal, iCell) = gradSSHZonal(iCell) SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) @@ -1031,7 +1033,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupCur, tracersGroupProvis, tracersGroupTend - integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop logical, pointer :: config_use_tracerGroup type (mpas_pool_iterator_type) :: groupItr @@ -1071,13 +1073,15 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) call mpas_pool_get_array(tendPool, 'lowFreqDivergence', lowFreqDivergenceTend) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = 1, maxLevelEdgeTop(iEdge) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) normalVelocityProvis(k, iEdge) = normalVelocityCur(k, iEdge) + rkSubstepWeight & * normalVelocityTend(k, iEdge) normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge)) @@ -1087,7 +1091,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessProvis(k, iCell) = layerThicknessCur(k, iCell) + rkSubstepWeight & * layerThicknessTend(k, iCell) end do @@ -1111,7 +1115,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersGroupProvis(:, k, iCell) = ( layerThicknessCur(k, iCell) * tracersGroupCur(:, k, iCell) & + rkSubstepWeight * tracersGroupTend(:, k, iCell) & ) / layerThicknessProvis(k, iCell) @@ -1208,7 +1212,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew, tracersGroupTend - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell logical, pointer :: config_use_tracerGroup type (mpas_pool_iterator_type) :: groupItr @@ -1243,12 +1247,13 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ call mpas_pool_get_array(tendPool, 'highFreqThickness', highFreqThicknessTend) call mpas_pool_get_array(tendPool, 'lowFreqDivergence', lowFreqDivergenceTend) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) layerThicknessNew(k, iCell) = layerThicknessNew(k, iCell) + rkWeight * layerThicknessTend(k, iCell) end do end do @@ -1277,7 +1282,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersGroupNew(:, k, iCell) = tracersGroupNew(:, k, iCell) + rkWeight & * tracersGroupTend(:, k, iCell) end do @@ -1326,7 +1331,7 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ real (kind=RKIND), dimension(:, :), pointer :: layerThicknessNew, normalVelocityNew real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew - integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: minLevelCell, maxLevelCell logical, pointer :: config_use_tracerGroup type (mpas_pool_iterator_type) :: groupItr @@ -1355,6 +1360,7 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_begin_iteration(tracersPool) @@ -1365,7 +1371,7 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) + do k = minLevelCell(iCell), maxLevelCell(iCell) tracersGroupNew(:, k, iCell) = tracersGroupNew(:, k, iCell) / layerThicknessNew(k, iCell) end do end do From 74d2acb836b79e20f3676993c4b937b813386fbe Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Mar 2021 08:54:07 -0600 Subject: [PATCH 41/46] Revert minLevel* changes to RK4: initialize *new = *cur --- src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F index f4e2d6f676..ce033f908d 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F @@ -242,7 +242,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + do k = 1, maxLevelEdgeTop(iEdge) normalVelocityNew(k, iEdge) = normalVelocityCur(k, iEdge) end do end do @@ -250,7 +250,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = minLevelCell(iCell), maxLevelCell(iCell) + do k = 1, maxLevelCell(iCell) layerThicknessNew(k, iCell) = layerThicknessCur(k, iCell) end do end do @@ -269,7 +269,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells ! couple tracers to thickness - do k = minLevelCell(iCell), maxLevelCell(iCell) + do k = 1, maxLevelCell(iCell) tracersNew(:, k, iCell) = tracersCur(:, k, iCell) * layerThicknessCur(k, iCell) end do end do From 6da36d2afe1e0a592637f2a8ea5014eed43b7f2d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Mar 2021 09:17:19 -0600 Subject: [PATCH 42/46] Revert minLevel* changes to RK4: diagnostic update --- src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F index ce033f908d..2581f6cbac 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F @@ -1081,7 +1081,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges - do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + do k = 1, maxLevelEdgeTop(iEdge) normalVelocityProvis(k, iEdge) = normalVelocityCur(k, iEdge) + rkSubstepWeight & * normalVelocityTend(k, iEdge) normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge)) @@ -1091,7 +1091,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = minLevelCell(iCell), maxLevelCell(iCell) + do k = 1, maxLevelCell(iCell) layerThicknessProvis(k, iCell) = layerThicknessCur(k, iCell) + rkSubstepWeight & * layerThicknessTend(k, iCell) end do @@ -1115,7 +1115,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells - do k = minLevelCell(iCell), maxLevelCell(iCell) + do k = 1, maxLevelCell(iCell) tracersGroupProvis(:, k, iCell) = ( layerThicknessCur(k, iCell) * tracersGroupCur(:, k, iCell) & + rkSubstepWeight * tracersGroupTend(:, k, iCell) & ) / layerThicknessProvis(k, iCell) From f0c76b7d7f38d8a657d8a8dc09c14ad9e99b2909 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Apr 2021 10:32:23 -0600 Subject: [PATCH 43/46] Fixup RK4: Add subpool calls Revert RK4 New=Cur subpool retrievals for minLevel* Use minLevelCell for velocityZonal,Meridional --- .../mode_forward/mpas_ocn_time_integration_rk4.F | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F index 2581f6cbac..9d98aa1e49 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F @@ -141,7 +141,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ ! Diagnostics Indices ! Mesh array pointers - integer, dimension(:), pointer :: minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelCell, maxLevelCell, maxLevelEdgeTop real (kind=RKIND), dimension(:), pointer :: bottomDepth ! Provis Array Pointers @@ -234,11 +234,12 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergenceCur, 1) call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergenceNew, 2) - call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + ! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur + ! assignments below are needed to maintain bit-for-bit results + !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges @@ -622,6 +623,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityNew, 2) @@ -721,8 +724,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells - surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(minLevelEdgeBot(iEdge), iCell) - surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(minLevelEdgeBot(iEdge), iCell) + surfaceVelocity(indexSurfaceVelocityZonal, iCell) = velocityZonal(minLevelCell(iCell), iCell) + surfaceVelocity(indexSurfaceVelocityMeridional, iCell) = velocityMeridional(minLevelCell(iCell), iCell) SSHGradient(indexSSHGradientZonal, iCell) = gradSSHZonal(iCell) SSHGradient(indexSSHGradientMeridional, iCell) = gradSSHMeridional(iCell) @@ -1033,7 +1036,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupCur, tracersGroupProvis, tracersGroupTend - integer, dimension(:), pointer :: minLevelCell, maxLevelCell, minLevelEdgeBot, maxLevelEdgeTop + integer, dimension(:), pointer :: minLevelCell, maxLevelCell, maxLevelEdgeTop logical, pointer :: config_use_tracerGroup type (mpas_pool_iterator_type) :: groupItr @@ -1075,7 +1078,6 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) !$omp parallel From 081693555f20a8b199b7522d151877029d94d3cd Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 9 Apr 2021 08:36:23 -0600 Subject: [PATCH 44/46] Keep same vertical limits on kpp call --- src/core_ocean/shared/mpas_ocn_vmix_cvmix.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F index bccd588920..4a6b064a94 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix_cvmix.F @@ -662,8 +662,8 @@ subroutine ocn_vmix_coefs_cvmix_build(meshPool, statePool, forcingPool, err, tim Snonlocal = cvmix_variables%kpp_Snonlocal_iface(minLevelCell(iCell):maxLevelCell(iCell)+1), & surf_fric = cvmix_variables%SurfaceFriction, & surf_buoy = cvmix_variables%SurfaceBuoyancyForcing, & - nlev = maxLevelCell(iCell) - minLevelCell(iCell) + 2, & - max_nlev = maxLevelCell(iCell) - minLevelCell(iCell) + 2, & + nlev = maxLevelCell(iCell) - minLevelCell(iCell) + 1, & + max_nlev = nVertLevels, & Langmuir_EFactor = langmuirEnhancementFactor ) ! call mpas_timer_stop('cvmix coeffs kpp') From 19dc7d971424324cf82e81262c6e4f102aee9536 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 12 Apr 2021 11:10:19 -0500 Subject: [PATCH 45/46] Fixup redi: change k to minLevelEdgeBot for flux terms --- src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F index eeb8a12faf..89b98f788e 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_hmix_redi.F @@ -248,13 +248,13 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg /(zMid(k, cell2) - zMid(k + 1, cell2))) if (minLevelCell(cell1) < minLevelEdgeBot(iEdge)) then flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & - layerThicknessEdge(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 1, iEdge)* & - (tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1))/(zMid(k - 1, cell1) - zMid(k, cell1)) + layerThicknessEdge(minLevelEdgeBot(iEdge), iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(minLevelEdgeBot(iEdge), 1, iEdge)* & + (tracers(iTr, minLevelEdgeBot(iEdge) - 1, cell1) - tracers(iTr, minLevelEdgeBot(iEdge), cell1))/(zMid(minLevelEdgeBot(iEdge) - 1, cell1) - zMid(minLevelEdgeBot(iEdge), cell1)) endif if (minLevelCell(cell2) < minLevelEdgeBot(iEdge)) then flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & - layerThicknessEdge(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 2, iEdge)* & - (tracers(iTr, k - 1, cell2) - tracers(iTr, k, cell2))/(zMid(k - 1, cell2) - zMid(k, cell2)) + layerThicknessEdge(minLevelEdgeBot(iEdge), iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(minLevelEdgeBot(iEdge), 2, iEdge)* & + (tracers(iTr, minLevelEdgeBot(iEdge) - 1, cell2) - tracers(iTr, minLevelEdgeBot(iEdge), cell2))/(zMid(minLevelEdgeBot(iEdge) - 1, cell2) - zMid(minLevelEdgeBot(iEdge), cell2)) endif redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2*invAreaCell From 9d189a2c2a810ec0e50652a67a11a07305ad0756 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 14 Apr 2021 11:42:29 -0600 Subject: [PATCH 46/46] Fixup ocn_tracer_ecosys --- .../shared/mpas_ocn_tracer_ecosys.F | 25 +++++++++---------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F b/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F index 9c2c166077..9e66758726 100755 --- a/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_ecosys.F @@ -1016,7 +1016,7 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, type (mpas_pool_type), pointer :: ecosysAuxiliary, & ! additional forcing fields ecosysSeaIceCoupling - integer :: numColumns, column, iCell, iTracer, iLevelSurface + integer :: numColumns, column, iCell, iTracer, kSurface ! input flux components in forcing pool real (kind=RKIND), dimension(:), pointer :: & @@ -1173,7 +1173,7 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, !DWJ 08/05/2016: This loop needs OpenMP added to it. do iCell=1,nCellsSolve - iLevelSurface = minLevelCell(iCell) + kSurface = minLevelCell(iCell) ! BGC_forcing%surfacePressure(column) = atmosphericPressure(iCell)*PascalsToAtmospheres BGC_forcing%surfacePressure(column) = 1.0_RKIND BGC_forcing%iceFraction(column) = iceFraction(iCell) @@ -1186,17 +1186,16 @@ subroutine ocn_tracer_ecosys_surface_flux_compute(activeTracers, ecosysTracers, BGC_forcing%atmCO2_ALT_CO2(column) = atmosphericCO2_ALT_CO2(iCell) BGC_forcing%surface_pH(column) = PH_PREV(iCell) BGC_forcing%surface_pH_alt_co2(column) = PH_PREV_ALT_CO2(iCell) - BGC_forcing%surfaceDepth(column) = -1.0_RKIND*zMid(iLevelSurface,iCell) - BGC_forcing%SST(column) = activeTracers(indexTemperature,iLevelSurface,iCell) - BGC_forcing%SSS(column) = activeTracers(indexSalinity,iLevelSurface,iCell) - - ! the 1 indices may need to be changed to iLevelSurface - BGC_input%BGC_tracers(1,column,BGC_indices%dic_ind) = ecosysTracers(ecosysIndices%dic_ind,1,iCell) - BGC_input%BGC_tracers(1,column,BGC_indices%dic_alt_co2_ind) = ecosysTracers(ecosysIndices%dic_alt_co2_ind,1,iCell) - BGC_input%BGC_tracers(1,column,BGC_indices%alk_ind) = ecosysTracers(ecosysIndices%alk_ind,1,iCell) - BGC_input%BGC_tracers(1,column,BGC_indices%po4_ind) = ecosysTracers(ecosysIndices%po4_ind,1,iCell) - BGC_input%BGC_tracers(1,column,BGC_indices%sio3_ind) = ecosysTracers(ecosysIndices%sio3_ind,1,iCell) - BGC_input%BGC_tracers(1,column,BGC_indices%o2_ind) = ecosysTracers(ecosysIndices%o2_ind,1,iCell) + BGC_forcing%surfaceDepth(column) = -1.0_RKIND*zMid(kSurface,iCell) + BGC_forcing%SST(column) = activeTracers(indexTemperature,kSurface,iCell) + BGC_forcing%SSS(column) = activeTracers(indexSalinity,kSurface,iCell) + + BGC_input%BGC_tracers(kSurface,column,BGC_indices%dic_ind) = ecosysTracers(ecosysIndices%dic_ind,kSurface,iCell) + BGC_input%BGC_tracers(kSurface,column,BGC_indices%dic_alt_co2_ind) = ecosysTracers(ecosysIndices%dic_alt_co2_ind,kSurface,iCell) + BGC_input%BGC_tracers(kSurface,column,BGC_indices%alk_ind) = ecosysTracers(ecosysIndices%alk_ind,kSurface,iCell) + BGC_input%BGC_tracers(kSurface,column,BGC_indices%po4_ind) = ecosysTracers(ecosysIndices%po4_ind,kSurface,iCell) + BGC_input%BGC_tracers(kSurface,column,BGC_indices%sio3_ind) = ecosysTracers(ecosysIndices%sio3_ind,kSurface,iCell) + BGC_input%BGC_tracers(kSurface,column,BGC_indices%o2_ind) = ecosysTracers(ecosysIndices%o2_ind,kSurface,iCell) ! NOTE pass in total Fe and mult by parm_Fe_bioavail inside the flux routine ! divide river Fe by bioavail since it is already the available to make it total