From 0723b5d99661bbfda6b82a945cb7b0060c04c2ee Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Wed, 27 Jul 2022 13:59:37 -0500 Subject: [PATCH 01/28] Adds submesoscale parameterization Adds the Fox-Kemper et al 2011 parameterization for submesocale eddies. There is also a substantial eddy parameterization reorganization as the buoyancy gradient that was calculated within GM is also needed for this parameterization. There are also cases when we will run with the submeso parameterization and not GM so that calculation is separate. The density threshold MLD is also removed from the analysis member as it is required by numerous calculations in the forward model now. The contributions of submesoscale eddy velocity to AMOC is also added to the MOC streamfunction computation --- components/mpas-ocean/bld/build-namelist | 21 +- .../mpas-ocean/bld/build-namelist-group-list | 2 + .../mpas-ocean/bld/build-namelist-section | 19 +- .../namelist_defaults_mpaso.xml | 13 +- .../namelist_definition_mpaso.xml | 78 ++++- .../mpas-ocean/cime_config/config_pes.xml | 36 +-- components/mpas-ocean/driver/ocn_comp_mct.F | 12 + components/mpas-ocean/src/Registry.xml | 72 ++++- .../Registry_mixed_layer_depths.xml | 12 - .../mpas_ocn_mixed_layer_depths.F | 73 +---- .../mpas_ocn_mixed_layer_heat_budget.F | 13 +- .../mpas_ocn_moc_streamfunction.F | 27 +- .../src/mode_forward/mpas_ocn_forward_mode.F | 12 + .../mpas_ocn_time_integration_rk4.F | 22 +- .../mpas_ocn_time_integration_si.F | 12 +- .../mpas_ocn_time_integration_split.F | 13 +- components/mpas-ocean/src/ocean.cmake | 2 + components/mpas-ocean/src/shared/Makefile | 12 +- .../src/shared/mpas_ocn_diagnostics.F | 38 ++- .../shared/mpas_ocn_diagnostics_variables.F | 78 ++++- .../mpas_ocn_eddy_parameterization_helpers.F | 282 ++++++++++++++++++ .../mpas-ocean/src/shared/mpas_ocn_gm.F | 131 +------- .../src/shared/mpas_ocn_submesoscale_eddies.F | 181 +++++++++++ .../src/shared/mpas_ocn_tracer_hmix_redi.F | 26 +- 24 files changed, 874 insertions(+), 313 deletions(-) create mode 100644 components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F create mode 100644 components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 64a922e9781d..6011661e6ef3 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -568,6 +568,16 @@ add_default($nl, 'config_Redi_horizontal_taper'); add_default($nl, 'config_Redi_horizontal_ramp_min'); add_default($nl, 'config_Redi_horizontal_ramp_max'); +###################################################### +# Namelist group: submesoscale_eddy_parameterization # +###################################################### + +add_default($nl, 'config_submesoscale_enable'); +add_default($nl, 'config_submesoscale_tau'); +add_default($nl, 'config_submesoscale_Ce'); +add_default($nl, 'config_submesoscale_Lfmin'); +add_default($nl, 'config_submesoscale_ds_max'); + ############################################ # Namelist group: GM_eddy_parameterization # ############################################ @@ -592,6 +602,13 @@ add_default($nl, 'config_GM_horizontal_ramp_max'); add_default($nl, 'config_GMRedi_Rossby_ramp_min'); add_default($nl, 'config_GMRedi_Rossby_ramp_max'); +######################################### +# Namelist group: eddy_parameterization # +######################################### + +add_default($nl, 'config_mixedLayerDepths_crit_dens_threshold'); +add_default($nl, 'config_mld_reference_depth'); + #################################### # Namelist group: Rayleigh_damping # #################################### @@ -1291,9 +1308,7 @@ add_default($nl, 'config_AM_mixedLayerDepths_output_stream'); add_default($nl, 'config_AM_mixedLayerDepths_write_on_startup'); add_default($nl, 'config_AM_mixedLayerDepths_compute_on_startup'); add_default($nl, 'config_AM_mixedLayerDepths_Tthreshold'); -add_default($nl, 'config_AM_mixedLayerDepths_Dthreshold'); add_default($nl, 'config_AM_mixedLayerDepths_crit_temp_threshold'); -add_default($nl, 'config_AM_mixedLayerDepths_crit_dens_threshold'); add_default($nl, 'config_AM_mixedLayerDepths_reference_pressure'); add_default($nl, 'config_AM_mixedLayerDepths_Tgradient'); add_default($nl, 'config_AM_mixedLayerDepths_Dgradient'); @@ -1681,7 +1696,9 @@ my @groups = qw(run_modes hmix_del4 hmix_leith redi_isopycnal_mixing + submesoscale_eddy_parameterization gm_eddy_parameterization + eddy_parameterization rayleigh_damping cvmix gotm diff --git a/components/mpas-ocean/bld/build-namelist-group-list b/components/mpas-ocean/bld/build-namelist-group-list index 6c332ff45f01..d71815ebbcb7 100644 --- a/components/mpas-ocean/bld/build-namelist-group-list +++ b/components/mpas-ocean/bld/build-namelist-group-list @@ -9,7 +9,9 @@ my @groups = qw(run_modes hmix_del4 hmix_leith redi_isopycnal_mixing + submesoscale_eddy_parameterization gm_eddy_parameterization + eddy_parameterization rayleigh_damping cvmix gotm diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index 0a9d5459111c..c1c185840e24 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -111,6 +111,16 @@ add_default($nl, 'config_Redi_horizontal_taper'); add_default($nl, 'config_Redi_horizontal_ramp_min'); add_default($nl, 'config_Redi_horizontal_ramp_max'); +###################################################### +# Namelist group: submesoscale_eddy_parameterization # +###################################################### + +add_default($nl, 'config_submesoscale_enable'); +add_default($nl, 'config_submesoscale_tau'); +add_default($nl, 'config_submesoscale_Ce'); +add_default($nl, 'config_submesoscale_Lfmin'); +add_default($nl, 'config_submesoscale_ds_max'); + ############################################ # Namelist group: GM_eddy_parameterization # ############################################ @@ -135,6 +145,13 @@ add_default($nl, 'config_GM_horizontal_ramp_max'); add_default($nl, 'config_GMRedi_Rossby_ramp_min'); add_default($nl, 'config_GMRedi_Rossby_ramp_max'); +######################################### +# Namelist group: eddy_parameterization # +######################################### + +add_default($nl, 'config_mixedLayerDepths_crit_dens_threshold'); +add_default($nl, 'config_mld_reference_depth'); + #################################### # Namelist group: Rayleigh_damping # #################################### @@ -763,9 +780,7 @@ add_default($nl, 'config_AM_mixedLayerDepths_output_stream'); add_default($nl, 'config_AM_mixedLayerDepths_write_on_startup'); add_default($nl, 'config_AM_mixedLayerDepths_compute_on_startup'); add_default($nl, 'config_AM_mixedLayerDepths_Tthreshold'); -add_default($nl, 'config_AM_mixedLayerDepths_Dthreshold'); add_default($nl, 'config_AM_mixedLayerDepths_crit_temp_threshold'); -add_default($nl, 'config_AM_mixedLayerDepths_crit_dens_threshold'); add_default($nl, 'config_AM_mixedLayerDepths_reference_pressure'); add_default($nl, 'config_AM_mixedLayerDepths_Tgradient'); add_default($nl, 'config_AM_mixedLayerDepths_Dgradient'); diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index a1c31bf79ab2..35679ddd12f1 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -152,6 +152,13 @@ 30e3 40e3 + +.true. +172800 +0.06 +1000.0 +100000.0 + .true. .false. @@ -191,6 +198,10 @@ 0.5 3.0 + +0.03 +10 + .false. 0.0 @@ -725,9 +736,7 @@ .false. .true. .true. -.true. 0.2 -0.03 1.0E5 .false. .false. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 841de7c4e7fc..3943d0e501ae 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -499,6 +499,49 @@ Default: Defined in namelist_defaults.xml + + + +flag to enable the FK2011 parameterization for submesoscale eddies + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +timescale for frictional slumping of front (in seconds) + +Valid values: positive reals, between 1-10 days +Default: Defined in namelist_defaults.xml + + + +efficiency of submesoscale eddies + +Valid values: 0.06 - 0.08 +Default: Defined in namelist_defaults.xml + + + +minimum frontal width (meters) + +Valid values: between 500 and 2000m +Default: Defined in namelist_defaults.xml + + + +maximum grid scale to scale up buoyancy gradient + +Valid values: around 1 degree +Default: Defined in namelist_defaults.xml + + + + + + +potential density change relative to surface for mixed layer depth threshold method. This calculation is used for the Redi tapering, GM N2_dependent bolus kappa, and the submesoscale eddy parameterization + +Valid values: suggested range 0.01 less than or equal to thresh less than or equal to 0.5 +Default: Defined in namelist_defaults.xml + + + +reference depth for threshold computation + +Valid values: any positive real, near 10 +Default: Defined in namelist_defaults.xml + + + - -Logical flag that determines if MLDs are calculated using a critical density threshold - -Valid values: .true. or .false. -Default: Defined in namelist_defaults.xml - - temperature change relative to surface for threshold method @@ -4036,14 +4090,6 @@ Valid values: all real positive values, suggested range 0.2 less than or equal t Default: Defined in namelist_defaults.xml - -potential density change relative to surface for threshold method - -Valid values: suggested range 0.01 less than or equal to thresh less than or equal to 0.5 -Default: Defined in namelist_defaults.xml - - reference pressure for threshold computation diff --git a/components/mpas-ocean/cime_config/config_pes.xml b/components/mpas-ocean/cime_config/config_pes.xml index 7cd87daed67d..dcdf96722438 100644 --- a/components/mpas-ocean/cime_config/config_pes.xml +++ b/components/mpas-ocean/cime_config/config_pes.xml @@ -235,34 +235,34 @@ mpas-ocean+chrysalis: standard-res, compset=DATM+MPASO, 15x32x2 NODESxMPIxOMP - 32 + 64 64 - 160 - 160 - 160 - 160 - 320 - 160 - 160 - 160 + 320 + 320 + 320 + 320 + 640 + 320 + 320 + 320 - 2 - 2 - 2 - 2 - 2 - 2 - 2 - 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 0 0 0 0 - 160 + 320 0 0 0 diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index 0156c31e1dcb..3517de9bafbb 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -60,6 +60,8 @@ module ocn_comp_mct use ocn_tracer_surface_restoring use ocn_gm use ocn_config + use ocn_submesoscale_eddies + use ocn_eddy_parameterization_helpers ! ! !PUBLIC MEMBER FUNCTIONS: implicit none @@ -971,6 +973,16 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr) + call ocn_eddy_compute_mixed_layer_depth(statePool) + if (config_use_GM.or.config_submesoscale_enable) then + call ocn_eddy_compute_buoyancy_gradient() + end if + + if (config_submesoscale_enable) then + call mpas_timer_start("submesoscale eddy velocity compute", .false.) + call ocn_submesoscale_compute_velocity() + call mpas_timer_stop("submesoscale eddy velocity compute") + end if ! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on if (config_use_Redi.or.config_use_GM) then call ocn_gm_compute_Bolus_velocity(statePool, meshPool, scratchPool, timeLevelIn=1) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index f8dc85f72e1f..1515ad1b2f59 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -343,6 +343,28 @@ possible_values="Any positive real value." /> + + + + + + + + + + + + + - + + + + + + @@ -2922,6 +2986,12 @@ + + - - - @@ -101,7 +90,6 @@ - diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_depths.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_depths.F index 7c3a85962104..ae951f533915 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_depths.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_depths.F @@ -174,13 +174,13 @@ subroutine ocn_compute_mixed_layer_depths(domain, timeLevel, err)!{{{ real (kind=RKIND), dimension(:,:), allocatable :: pressureAdjustedForLandIce real (kind=RKIND), dimension(:), pointer :: tThreshMLD, tGradientMLD, landIceDraft - real (kind=RKIND), dimension(:), pointer :: dThreshMLD, dGradientMLD, landIcePressure + real (kind=RKIND), dimension(:), pointer :: dGradientMLD, landIcePressure integer, dimension(:), pointer :: landIceMask real (kind=RKIND), dimension(:,:,:), pointer :: tracers integer :: interp_local real (kind=RKIND), allocatable, dimension(:,:) :: densityGradient, temperatureGradient - real (kind=RKIND) :: mldTemp,dTempThres, dDenThres, dTempGrad, dDenGrad - real (kind=RKIND) :: dz,temp_ref_lev, den_ref_lev, dV, dVm1, dVp1, localVals(6) + real (kind=RKIND) :: mldTemp,dTempThres, dTempGrad, dDenGrad + real (kind=RKIND) :: dz,temp_ref_lev, dV, dVm1, dVp1, localVals(6) real (kind=RKIND), dimension(:), pointer :: latCell, lonCell err = 0 @@ -294,73 +294,6 @@ subroutine ocn_compute_mixed_layer_depths(domain, timeLevel, err)!{{{ end if !end tThresh MLD search - if(config_AM_mixedLayerDepths_Dthreshold) then - call mpas_pool_get_array(mixedLayerDepthsAMPool, 'dThreshMLD',dThreshMLD) - - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(refIndex, found_den_mld, k, localvals, den_ref_lev, dVp1, dV, mldTemp) - do iCell = 1,nCells - - !Initialize RefIndex for cases of very shallow columns - refIndex = maxLevelCell(iCell) - found_den_mld = .false. - - do k=1, maxLevelCell(iCell)-1 - if(pressureAdjustedForLandIce(k+1,iCell) > config_AM_mixedLayerDepths_reference_pressure) then - localvals(2:3)=potentialDensity(k:k+1,iCell) - localvals(5:6)=pressureAdjustedForLandIce(k:k+1,iCell) - - call interp_bw_levels(localVals(2),localVals(3),localVals(5),localVals(6), & - config_AM_mixedLayerDepths_reference_pressure,interp_local,den_ref_lev) - refIndex=k - exit - end if - end do - - do k=refIndex,maxLevelCell(iCell)-1 - - if(.not. found_den_mld .and. abs(potentialDensity(k+1,iCell) - den_ref_lev) .ge. & - config_AM_mixedLayerDepths_crit_dens_threshold) then - dVp1 = abs(potentialDensity(k+1,iCell) - den_ref_lev) - dV = abs(potentialDensity(k ,iCell) - den_ref_lev) - localVals(2:3)=zMid(k:k+1,iCell) - call interp_bw_levels(localVals(2),localVals(3), dV, dVp1, config_AM_mixedLayerDepths_crit_dens_threshold, & - interp_local, mldTemp)!, dVm1, localVals(1)) - mldTemp=max(mldTemp,zMid(k+1,iCell)) !make sure MLD isn't deeper than zMid(k+1) - dThreshMLD(iCell)=abs(min(mldTemp,zMid(k,iCell))) !MLD should be deeper than zMid(k) - indMLD(iCell) = k+1 - found_den_mld = .true. - exit - end if - end do - -! if the mixed layer depth is not found, it is set to the depth of the bottom most level - - if(.not. found_den_mld) then - dThreshMLD(iCell) = abs(zMid(maxLevelCell(iCell),iCell)) - indMLD(iCell) = maxLevelCell(iCell) - end if - end do !iCell - !$omp end do - !$omp end parallel - - if (associated(landIceMask) ) then - !$omp parallel - !$omp do schedule(runtime) - do iCell = 1, nCells - dThreshMLD(iCell) = dThreshMLD(iCell) - abs(landIceDraft(iCell))*landIceMask(iCell) - end do - !$omp end do - !$omp end parallel - - - end if - - - end if !end dThresh MLD search - - ! Compute the mixed layer depth based on a gradient threshold in temperature and density if(config_AM_mixedLayerDepths_Tgradient) then call mpas_pool_get_array(mixedLayerDepthsAMPool, 'tGradMLD',tGradientMLD) diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_heat_budget.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_heat_budget.F index f27b593833b9..07648b39be70 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_heat_budget.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_mixed_layer_heat_budget.F @@ -103,18 +103,8 @@ subroutine ocn_init_mixed_layer_heat_budget(domain, err) !{{{ ! !----------------------------------------------------------------- - logical, pointer :: dThresholdFlag, mldON err = 0 - call mpas_pool_get_config(domain%configs, 'config_AM_mixedLayerDepths_Dthreshold', dThresholdFlag) - call mpas_pool_get_config(domain%configs, 'config_AM_mixedLayerDepths_enable', mldON) - - if (.not. mldON .and. .not. dThresholdFlag) THEN - call mpas_log_write("ML Heat Budget AM requires config_AM_mixedLayerDepths_enable and " & - //"config_AM_mixedLayerDepths_Dthreshold to both be true", MPAS_LOG_CRIT) - err = 1 - endif - end subroutine ocn_init_mixed_layer_heat_budget !}}} !*********************************************************************** @@ -184,7 +174,7 @@ subroutine ocn_compute_mixed_layer_heat_budget(domain, timeLevel, err) !{{{ activeTracerNonLocalMLTend, activeTracerHorMixMLTend, activeTracerVertMixMLTend, & activeTracerHorAdvectionMLTend, activeTracerVertAdvectionMLTend, activeTracersML, & layerThickness, activeTracersTendML - real(kind=RKIND), dimension(:), pointer :: dThreshMLD, bruntVaisalaFreqML + real(kind=RKIND), dimension(:), pointer :: bruntVaisalaFreqML real(kind=RKIND) :: depth, difference err = 0 @@ -222,7 +212,6 @@ subroutine ocn_compute_mixed_layer_heat_budget(domain, timeLevel, err) !{{{ call mpas_pool_get_array(mixedLayerHeatBudgetAMPool, 'activeTracersML', activeTracersML) call mpas_pool_get_array(mixedLayerHeatBudgetAMPool, 'activeTracersTendML', activeTracersTendML) call mpas_pool_get_array(mixedLayerHeatBudgetAMPool, 'bruntVaisalaFreqML', bruntVaisalaFreqML) - call mpas_pool_get_array(mixedLayerDepthsAMPool, 'dThreshMLD', dThreshMLD) call mpas_pool_get_dimension(tracersPool, 'index_temperature', index_temperature) call mpas_pool_get_array(tracersPool, 'activeTracers', tracers, timeLevel) call mpas_pool_get_array(tracerTendPool, 'activeTracersTend', tracerTend) diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F index 66dcc068d7f5..5a2f71d71995 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F @@ -507,10 +507,12 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ * 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2)) end do end do - do k = 2, nVertLevels + mocStreamValLatAndDepthRegionLocal(1, nVertLevels, iTransect) = & + - sumTransport(nVertLevels, iTransect) + do k = nVertLevels-1,1,-1 mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & - mocStreamValLatAndDepthRegionLocal(1, k - 1, iTransect) & - + sumTransport(k - 1, iTransect) + mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & + - sumTransport(k + 1, iTransect) end do end do @@ -602,14 +604,18 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ sumTransport(k,iTransect) = sumTransport(k,iTransect) + & transectEdgeMaskSigns(currentTransect,iEdge) & * transectEdgeMasks(currentTransect, iEdge) & - * normalVelocity(k,iEdge)*dvEdge(iEdge) & + * (normalMLEvelocity(k,iEdge) + normalVelocity(k,iEdge)) & + *dvEdge(iEdge) & * 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2)) end do end do - do k = 2, nVertLevels + ! do k = 2, nVertLevels + + mocStreamValLatAndDepthRegionLocal(1, nVertLevels, iTransect) = -sumTransport(nVertLevels,iTransect) + do k = nVertLevels-1,1,-1 mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & - mocStreamValLatAndDepthRegionLocal(1, k - 1, iTransect) & - + sumTransport(k - 1, iTransect) + mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & + - sumTransport(k + 1, iTransect) end do end do @@ -620,10 +626,13 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ do k = 1,maxLevelCell(iCell) do i = 1, regionsInAddGroup currentRegion = regionsInGroup(i, regionGroupNumber) - sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + (vertGMBolusVelocityTop(k, iCell) * & + sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + & + ((vertGMBolusVelocityTop(k, iCell) +vertMLEBolusVelocityTop(k,iCell)) * & areaCell(iCell) * regionCellMasks(currentRegion, iCell)) end do - sumVertBinVelocity(iBin, k) = sumVertBinVelocity(iBin, k) + (vertGMBolusVelocityTop(k, iCell) * areaCell(iCell)) + sumVertBinVelocity(iBin, k) = sumVertBinVelocity(iBin, k) + & + ((vertGMBolusVelocityTop(k, iCell) + vertMLEBolusVelocityTop(k,iCell)) & + * areaCell(iCell)) end do end do diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F index d007eaedd8c3..f8d69d6ff380 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F @@ -43,6 +43,7 @@ module ocn_forward_mode use ocn_test use ocn_mesh + use ocn_eddy_parameterization_helpers use ocn_thick_hadv use ocn_thick_vadv use ocn_thick_ale @@ -75,6 +76,7 @@ module ocn_forward_mode use ocn_tracer_CFC use ocn_tracer_surface_restoring use ocn_gm + use ocn_submesoscale_eddies use ocn_high_freq_thickness_hmix_del2 @@ -441,6 +443,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ ierr = ior(ierr,err_tmp) call ocn_gm_init(domain, err_tmp) ierr = ior(ierr,err_tmp) + call ocn_submesoscale_init(err_tmp) + ierr = ior(ierr,err_tmp) call ocn_tracer_nonlocalflux_init(err_tmp) ierr = ior(ierr,err_tmp) call ocn_tracer_ecosys_init(domain, err_tmp) @@ -675,6 +679,14 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{ call ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePool, err) ! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on + + call ocn_eddy_compute_mixed_layer_depth(statePool) + if (config_use_GM.or.config_submesoscale_enable) then + call mpas_timer_start("submesoscale eddy velocity compute", .false.) + call ocn_eddy_compute_buoyancy_gradient() + call mpas_timer_stop("submesoscale eddy velocity compute") + end if + if (config_use_Redi.or.config_use_GM) then call ocn_gm_compute_Bolus_velocity(statePool, meshPool, scratchPool, timeLevelIn=1) end if diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 844b5bc2483b..454ae89ac1eb 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -711,7 +711,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessNew, normalVelocityNew) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity, normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) @@ -784,7 +784,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ !$omp parallel !$omp do schedule(runtime) do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + normalGMBolusVelocity(:, iEdge) + normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + & + normalGMBolusVelocity(:, iEdge) + & + normalMLEvelocity(:, iEdge) end do !$omp end do !$omp end parallel @@ -950,7 +952,9 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh !$omp parallel !$omp do schedule(runtime) do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + normalGMBolusVelocity(:, iEdge) + normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + & + normalGMBolusVelocity(:, iEdge) + & + normalMLEvelocity(:, iEdge) end do !$omp end do !$omp end parallel @@ -1020,7 +1024,9 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig !$omp parallel !$omp do schedule(runtime) do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + normalGMBolusVelocity(:,iEdge) + normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + & + normalGMBolusVelocity(:,iEdge) + & + normalMLEvelocity(:,iEdge) end do !$omp end do !$omp end parallel @@ -1299,7 +1305,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessCur, normalVelocityCur) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity,normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshCur) !$acc enter data copyin(activeTracersCur) @@ -1368,7 +1374,9 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, !$omp parallel !$omp do schedule(runtime) do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + normalGMBolusVelocity(:,iEdge) + normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + & + normalGMBolusVelocity(:,iEdge) + & + normalMLEvelocity(:,iEdge) end do !$omp end do !$omp end parallel @@ -1589,7 +1597,7 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessNew, normalVelocityNew) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity,normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 0b900c04702c..cda0954fd1ef 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -340,7 +340,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! barotropicForcing, barotropicThicknessFlux ! layerThickEdgeFlux, layerThickEdgeMean, ! normalTransportVelocity, normalGMBolusVelocity - ! vertAleTransportTop + ! vertAleTransportTop, normalMLEvelocity ! velocityX, velocityY, velocityZ ! velocityZonal, velocityMeridional ! gradSSH, gradSSHX, gradSSHY, gradSSHZ @@ -1677,7 +1677,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! normalBaroclinicVelocity + uBolus uTemp(:) = normalBarotropicVelocityNew(iEdge) & + normalBaroclinicVelocityNew(:,iEdge) & - + normalGMBolusVelocity(:,iEdge) + + normalGMBolusVelocity(:,iEdge) & + + normalMLEvelocity(:,iEdge) ! thicknessSum is initialized outside the loop because ! on land boundaries maxLevelEdgeTop=0, but I want to @@ -1717,6 +1718,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ *(normalBarotropicVelocityNew(iEdge) + & normalBaroclinicVelocityNew(k,iEdge) + & normalGMBolusVelocity(k,iEdge) + & + normalMLEvelocity(k,iEdge) + & normalVelocityCorrection ) enddo @@ -1875,7 +1877,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessNew, normalVelocityNew) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity, normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) @@ -2276,7 +2278,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessNew, normalVelocityNew) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity,normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) @@ -2402,7 +2404,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$acc enter data copyin(layerThicknessCur) endif !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity,normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(tracersGroupNew) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index d89847e37bf6..556a73dbd6b2 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -275,7 +275,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! barotropicForcing, barotropicThicknessFlux ! layerThickEdgeFlux, layerThickEdgeMean, ! normalTransportVelocity, normalGMBolusVelocity - ! vertAleTransportTop + ! vertAleTransportTop, normalMLEvelocity ! velocityX, velocityY, velocityZ ! velocityZonal, velocityMeridional ! gradSSH, gradSSHX, gradSSHY, gradSSHZ @@ -837,7 +837,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! in tendency calls in stage 3. normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge) & *(normalBaroclinicVelocityNew(k,iEdge) & - + normalGMBolusVelocity(k,iEdge)) + + normalGMBolusVelocity(k,iEdge) & + + normalMLEvelocity(k,iEdge)) enddo ! vertical end do ! iEdge @@ -1562,7 +1563,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! normalBaroclinicVelocity + uBolus uTemp(:) = normalBarotropicVelocityNew(iEdge) & + normalBaroclinicVelocityNew(:,iEdge) & - + normalGMBolusVelocity(:,iEdge) + + normalGMBolusVelocity(:,iEdge) & + + normalMLEvelocity(:,iEdge) ! thicknessSum is initialized outside the loop because ! on land boundaries maxLevelEdgeTop=0, but I want to @@ -1600,6 +1602,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ *(normalBarotropicVelocityNew(iEdge) + & normalBaroclinicVelocityNew(k,iEdge) + & normalGMBolusVelocity(k,iEdge) + & + normalMLEvelocity(k,iEdge) + & normalVelocityCorrection ) enddo @@ -2322,7 +2325,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #ifdef MPAS_OPENACC !$acc enter data copyin(layerThicknessNew, normalVelocityNew) !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity, normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) @@ -2466,7 +2469,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$acc enter data copyin(layerThicknessCur) endif !$acc update device (normalTransportVelocity, & - !$acc normalGMBolusVelocity) + !$acc normalGMBolusVelocity, normalMLEvelocity) !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshNew) !$acc enter data copyin(activeTracersNew) diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index fb1dba45d4ba..2bd3c6d274f8 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -32,6 +32,8 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_init_routines.F core_ocean/shared/mpas_ocn_gm.F + core_ocean/shared/mpas_ocn_submesoscale_eddies.F + core_ocean/shared/mpas_ocn_eddy_parameterization_helpers.F core_ocean/shared/mpas_ocn_diagnostics.F core_ocean/shared/mpas_ocn_diagnostics_variables.F core_ocean/shared/mpas_ocn_mesh.F diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 81d8235bb707..a76476d2ef2d 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -2,6 +2,8 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_gm.o \ + mpas_ocn_eddy_parameterization_helpers.o \ + mpas_ocn_submesoscale_eddies.o \ mpas_ocn_diagnostics.o \ mpas_ocn_diagnostics_variables.o \ mpas_ocn_thick_ale.o \ @@ -76,11 +78,11 @@ OBJS = mpas_ocn_init_routines.o \ all: $(OBJS) -mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o +mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_submesoscale_eddies.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o -mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o +mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_submesoscale_eddies.o mpas_ocn_diagnostics_variables.o: mpas_ocn_config.o @@ -96,7 +98,11 @@ mpas_ocn_thick_vadv.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_thick_surface_flux.o: mpas_ocn_forcing.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o -mpas_ocn_gm.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o +mpas_ocn_gm.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_submesoscale_eddies.o + +mpas_ocn_submesoscale_eddies.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o + +mpas_ocn_eddy_parameterization_helpers.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 3117cf6314e6..1e79bd9d6cdc 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -232,13 +232,13 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo ! inputs: relativeVorticity, layerThickEdgeFlux, normalVelocity, normalTransportVelocity, ! normalGMBolusVelocity ! outputs: vertTransportVelocityTop, vertGMBolusVelocityTop, relativeVorticityCell, - ! divergence, kineticEnergyCell, tangentialVelocity + ! divergence, kineticEnergyCell, tangentialVelocity, vertMLEBolusVelocityTop ! in and out: vertVelocityTop call ocn_diagnostic_solve_vortVel(relativeVorticity, & layerThickEdgeFlux, normalVelocity, normalTransportVelocity, & - normalGMBolusVelocity, vertVelocityTop, vertTransportVelocityTop, & - vertGMBolusVelocityTop, relativeVorticityCell, divergence, & - kineticEnergyCell, tangentialVelocity) + normalGMBolusVelocity, normalMLEvelocity, vertVelocityTop, vertTransportVelocityTop, & + vertGMBolusVelocityTop, vertMLEBolusVelocityTop, relativeVorticityCell, & + divergence, kineticEnergyCell, tangentialVelocity) if ( configVertAdvMethod == vertAdvRemap ) then @@ -1563,9 +1563,9 @@ end subroutine ocn_diagnostic_solve_surfaceLayer !}}} !----------------------------------------------------------------------- subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & layerThicknessEdgeFlux, normalVelocity, normalTransportVelocity, & - normalGMBolusVelocity, vertVelocityTop, vertTransportVelocityTop, & - vertGMBolusVelocityTop, relativeVorticityCell, divergence, & - kineticEnergyCell, tangentialVelocity)!{{{ + normalGMBolusVelocity, normalMLEVelocity, vertVelocityTop, vertTransportVelocityTop, & + vertGMBolusVelocityTop, vertMLEBolusVelocityTop, relativeVorticityCell, & + divergence, kineticEnergyCell, tangentialVelocity)!{{{ implicit none @@ -1585,6 +1585,8 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & normalTransportVelocity real (kind=RKIND), dimension(:,:), intent(in) :: & normalGMBolusVelocity + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalMLEVelocity !----------------------------------------------------------------- ! @@ -1598,6 +1600,8 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & vertTransportVelocityTop real (kind=RKIND), dimension(:,:), intent(out) :: & vertGMBolusVelocityTop + real (kind=RKIND), dimension(:,:), intent(out) :: & + vertMLEBolusVelocityTop real (kind=RKIND), dimension(:,:), intent(out) :: & relativeVorticityCell real (kind=RKIND), dimension(:,:), intent(out) :: & @@ -1621,8 +1625,10 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp real (kind=RKIND), dimension(:), allocatable:: div_hu,div_huTransport,div_huGMBolus + real (kind=RKIND), dimension(:), allocatable:: div_huMLEBolus + #ifdef MPAS_OPENACC - !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus) + !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus, div_huMLEBolus) #endif ! Need owned cells for relativeVorticityCell @@ -1664,6 +1670,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & ! Compute divergence, kinetic energy, and vertical velocity ! allocate(div_hu(nVertLevels),div_huTransport(nVertLevels),div_huGMBolus(nVertLevels)) + allocate(div_huMLEBolus(nVertLevels)) #ifdef MPAS_OPENACC !$acc parallel loop gang vector & @@ -1672,7 +1679,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & !$acc maxLevelCell, normalVelocity, divergence, layerThicknessEdgeFlux, & !$acc normalTransportVelocity, normalGMBolusVelocity, vertVelocityTop, & !$acc vertTransportVelocityTop, vertGMBolusVelocityTop, invAreaCell, & - !$acc minLevelCell) & + !$acc minLevelCell,vertMLEBolusVelocityTop, normalMLEVelocity) & !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & !$acc dvEdge_temp, k, r_tmp) #else @@ -1687,6 +1694,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & div_hu(:) = 0.0_RKIND div_huTransport(:) = 0.0_RKIND div_huGMBolus(:) = 0.0_RKIND + div_huMLEBolus(:) = 0.0_RKIND invAreaCell1 = invAreaCell(iCell) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) @@ -1708,6 +1716,9 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & div_huGMBolus(k) = div_huGMBolus(k) & - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & * normalGMBolusVelocity(k, iEdge) * invAreaCell1 + div_huMLEBolus(k) = div_huMLEBolus(k) & + - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & + * normalMLEVelocity(k, iEdge) * invAreaCell1 end do end do ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. @@ -1717,6 +1728,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & vertVelocityTop(k,iCell) = vertVelocityTop(k+1,iCell) - div_hu(k) vertTransportVelocityTop(k,iCell) = vertTransportVelocityTop(k+1,iCell) - div_huTransport(k) vertGMBolusVelocityTop(k,iCell) = vertGMBolusVelocityTop(k+1,iCell) - div_huGMBolus(k) + vertMLEBolusVelocityTop(k,iCell) = vertMLEBolusVelocityTop(k+1,iCell) - div_huMLEBolus(k) end do end do #ifndef MPAS_OPENACC @@ -3115,6 +3127,14 @@ subroutine ocn_reconstruct_gm_vectors(meshPool) !{{{ GMBolusVelocityMeridional & ) + call mpas_reconstruct(meshPool, normalMLEvelocity, & + mleVelocityX, & + mleVelocityY, & + mleVelocityZ, & + mleVelocityZonal, & + mleVelocityMeridional & + ) + call mpas_reconstruct(meshPool, gmStreamFuncTopOfEdge, & GMStreamFuncX, & GMStreamFuncY, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index 3325960bda3e..644257970ff5 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -62,7 +62,9 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:), pointer :: vertVelocityTop real (kind=RKIND), dimension(:,:), pointer :: vertTransportVelocityTop real (kind=RKIND), dimension(:,:), pointer :: vertGMBolusVelocityTop + real (kind=RKIND), dimension(:,:), pointer :: vertMLEBolusVelocityTop real (kind=RKIND), dimension(:,:), pointer :: normalGMBolusVelocity + real (kind=RKIND), dimension(:,:), pointer :: normalMLEvelocity real (kind=RKIND), dimension(:,:), pointer :: normalTransportVelocity real (kind=RKIND), dimension(:,:), pointer :: RiTopOfCell real (kind=RKIND), dimension(:,:), pointer :: RiTopOfEdge @@ -99,7 +101,8 @@ module ocn_diagnostics_variables transportVelocityX, transportVelocityY, transportVelocityZ, transportVelocityZonal, & transportVelocityMeridional, GMBolusVelocityX, GMBolusVelocityY, & GMBolusVelocityZ, GMBolusVelocityZonal, GMBolusVelocityMeridional, gmStreamFuncTopOfEdge, & - GMStreamFuncX, GMStreamFuncY, GMStreamFuncZ, GMStreamFuncZonal, GMStreamFuncMeridional + GMStreamFuncX, GMStreamFuncY, GMStreamFuncZ, GMStreamFuncZonal, GMStreamFuncMeridional, & + mleVelocityX, mleVelocityY, mleVelocityZ, mleVelocityZonal, mleVelocityMeridional real (kind=RKIND), dimension(:,:), pointer :: rx1Edge real (kind=RKIND), dimension(:,:), pointer :: rx1Cell @@ -119,11 +122,13 @@ module ocn_diagnostics_variables real(kind=RKIND), dimension(:,:), pointer :: RediKappaScaling real(kind=RKIND), dimension(:,:), pointer :: RediKappaSfcTaper real(kind=RKIND), dimension(:,:), pointer :: k33 + real(kind=RKIND), dimension(:,:), pointer :: gradBuoyEddy real(kind=RKIND), dimension(:), pointer :: gmBolusKappa real(kind=RKIND), dimension(:), pointer :: cGMphaseSpeed real(kind=RKIND), dimension(:,:,:), pointer :: slopeTriadUp, slopeTriadDown - integer, dimension(:), pointer :: indMLD + integer, dimension(:), pointer :: indMLD, indMLDedge + real(kind=RKIND), dimension(:), pointer :: dThreshMLD real (kind=RKIND), dimension(:,:), pointer :: velocityX, velocityY, velocityZ real (kind=RKIND), dimension(:,:), pointer :: vertAleTransportTop @@ -265,6 +270,11 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e RossbyRadius) end if + if (config_use_GM .or. config_submesoscale_enable) then + call mpas_pool_get_array(diagnosticsPool, 'gradBuoyEddy', & + gradBuoyEddy) + end if + if ( config_compute_active_tracer_budgets ) then call mpas_pool_get_array(diagnosticsPool, & 'activeTracerSurfaceFluxTendency', & @@ -407,7 +417,27 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e GMStreamFuncZonal) call mpas_pool_get_array(diagnosticsPool, 'GMStreamFuncMeridional', & GMStreamFuncMeridional) - call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityX', & + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityX', & + mleVelocityX) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityY', & + mleVelocityY) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZ', & + mleVelocityZ) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZonal', & + mleVelocityZonal) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & + mleVelocityMeridional) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityX', & + mleVelocityX) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityY', & + mleVelocityY) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZ', & + mleVelocityZ) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZonal', & + mleVelocityZonal) + call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & + mleVelocityMeridional) + call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityX', & GMBolusVelocityX) call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityY', & GMBolusVelocityY) @@ -419,8 +449,12 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e GMBolusVelocityMeridional) call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', & normalGMBolusVelocity) + call mpas_pool_get_array(diagnosticsPool, 'normalMLEvelocity', & + normalMLEvelocity) call mpas_pool_get_array(diagnosticsPool, 'vertGMBolusVelocityTop', & vertGMBolusVelocityTop) + call mpas_pool_get_array(diagnosticsPool, 'vertMLEBolusVelocityTop', & + vertMLEBolusVelocityTop) call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdgeFlux', & layerThickEdgeFlux) @@ -589,6 +623,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e 'modifyLandIcePressureMask', & modifyLandIcePressureMask) call mpas_pool_get_array(diagnosticsPool, 'indMLD', indMLD) + call mpas_pool_get_array(diagnosticsPool, 'indMLDedge', indMLDedge) + call mpas_pool_get_array(diagnosticsPool, 'dThreshMLD', dThreshMLD) call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', & surfacePressure) @@ -688,6 +724,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & + !$acc mleVelocityX, & + !$acc mleVelocityY, & + !$acc mleVelocityZ, & !$acc GMBolusVelocityX, & !$acc GMBolusVelocityY, & !$acc GMBolusVelocityZ, & @@ -710,12 +749,15 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc BruntVaisalaFreqTop, & !$acc vertAleTransportTop, & !$acc GMBolusVelocityZonal, & + !$acc mleVelocityZonal, & !$acc bulkRichardsonNumber, & !$acc gmStreamFuncTopOfEdge, & !$acc indexSSHGradientZonal, & !$acc relativeVorticityCell, & !$acc normalGMBolusVelocity, & + !$acc normalMLEvelocity, & !$acc vertGMBolusVelocityTop, & + !$acc vertMLEBolusVelocityTop, & !$acc transportVelocityZonal, & !$acc GMStreamFuncMeridional, & !$acc indexSurfaceLayerDepth, & @@ -728,6 +770,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc bulkRichardsonNumberBuoy, & !$acc tracersSurfaceLayerValue, & !$acc GMBolusVelocityMeridional, & + !$acc mleVelocityMeridional, & !$acc bulkRichardsonNumberShear, & !$acc indexSurfaceVelocityZonal, & !$acc normalVelocitySurfaceLayer, & @@ -757,6 +800,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc RossbyRadius & !$acc ) end if + if (config_use_GM.or.config_submesoscale_enable) then + !$acc enter data create(gradBuoyEddy) + end if if ( config_compute_active_tracer_budgets ) then !$acc enter data create( & !$acc activeTracerHorMixTendency, & @@ -813,6 +859,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc zTop, & !$acc xtime, & !$acc indMLD, & + !$acc indMLDedge, & + !$acc dThreshMLD, & !$acc density, & !$acc betaEdge, & !$acc pressure, & @@ -928,6 +976,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & + !$acc mleVelocityX, & + !$acc mleVelocityY, & + !$acc mleVelocityZ, & !$acc GMBolusVelocityX, & !$acc GMBolusVelocityY, & !$acc GMBolusVelocityZ, & @@ -950,12 +1001,15 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc BruntVaisalaFreqTop, & !$acc vertAleTransportTop, & !$acc GMBolusVelocityZonal, & + !$acc mleVelocityZonal, & !$acc bulkRichardsonNumber, & !$acc gmStreamFuncTopOfEdge, & !$acc indexSSHGradientZonal, & !$acc relativeVorticityCell, & !$acc normalGMBolusVelocity, & + !$acc normalMLEvelocity, & !$acc vertGMBolusVelocityTop, & + !$acc vertMLEBolusVelocityTop, & !$acc transportVelocityZonal, & !$acc GMStreamFuncMeridional, & !$acc indexSurfaceLayerDepth, & @@ -967,6 +1021,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc vertTransportVelocityTop, & !$acc bulkRichardsonNumberBuoy, & !$acc tracersSurfaceLayerValue, & + !$acc mleVelocityMeridional, & !$acc GMBolusVelocityMeridional, & !$acc bulkRichardsonNumberShear, & !$acc indexSurfaceVelocityZonal, & @@ -997,6 +1052,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc RossbyRadius & !$acc ) end if + if (config_use_GM.or.config_submesoscale_enable) then + !$acc exit data delete(gradBuoyEddy) + end if if ( config_compute_active_tracer_budgets ) then !$acc exit data delete( & !$acc activeTracerHorMixTendency, & @@ -1053,6 +1111,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc zTop, & !$acc xtime, & !$acc indMLD, & + !$acc indMLDedge, & + !$acc dThreshMLD, & !$acc density, & !$acc betaEdge, & !$acc pressure, & @@ -1129,6 +1189,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ GMBolusVelocityX, & GMBolusVelocityY, & GMBolusVelocityZ, & + mleVelocityX, & + mleVelocityY, & + mleVelocityZ, & vertNonLocalFlux, & vertViscTopOfEdge, & gradSSHMeridional, & @@ -1148,12 +1211,15 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ BruntVaisalaFreqTop, & vertAleTransportTop, & GMBolusVelocityZonal, & + mleVelocityZonal, & bulkRichardsonNumber, & gmStreamFuncTopOfEdge, & indexSSHGradientZonal, & relativeVorticityCell, & normalGMBolusVelocity, & + normalMLEvelocity, & vertGMBolusVelocityTop, & + vertMLEBolusVelocityTop, & transportVelocityZonal, & GMStreamFuncMeridional, & indexSurfaceLayerDepth, & @@ -1165,6 +1231,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ vertTransportVelocityTop, & bulkRichardsonNumberBuoy, & tracersSurfaceLayerValue, & + mleVelocityMeridional, & GMBolusVelocityMeridional, & bulkRichardsonNumberShear, & indexSurfaceVelocityZonal, & @@ -1190,6 +1257,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ gmHorizontalTaper, & RossbyRadius) end if + if (config_use_GM.or.config_submesoscale_enable) then + nullify(gradBuoyEddy) + end if if ( config_compute_active_tracer_budgets ) then nullify(activeTracerHorMixTendency, & activeTracerVertMixTendency, & @@ -1242,6 +1312,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ zTop, & xtime, & indMLD, & + indMLDedge, & + dThreshMLD, & density, & betaEdge, & pressure, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F new file mode 100644 index 000000000000..5c1f12faedba --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F @@ -0,0 +1,282 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! + +module ocn_eddy_parameterization_helpers + + use mpas_pool_routines + use mpas_derived_types + use mpas_constants + use mpas_threading + use mpas_timer + + use ocn_constants + use ocn_config + use ocn_diagnostics_variables + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_eddy_compute_buoyancy_gradient, & + ocn_eddy_compute_mixed_layer_depth + + + !*********************************************************************** + + contains + + !*********************************************************************** + ! + ! routine ocn_eddy_compute_buoyancy_gradient + ! + !> \brief Computes fixed depth horizontal buoyancy gradient + !> \details + !> Computes the horizontal buoyancy gradient on fixed depth levels. + !> This quantity is extracted from mpas_ocn_gm as it is needed for both + !> the GM eddy parameterization and the submesoscale eddy parameterization + !*********************************************************************** + + subroutine ocn_eddy_compute_buoyancy_gradient() + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2, indMLDedge + ! gradDensityEdge: Normal gradient of density + ! units: none + real(kind=RKIND), dimension(:, :), allocatable :: gradDensityEdge, gradDensityTopOfEdge, & + dDensityDzTopOfCell, dDensityDzTopOfEdge, & + gradZMidEdge, gradZMidTopOfEdge + real(kind=RKIND) :: rtmp, h1, h2 + real(kind=RKIND), parameter :: epsGM = 1.0E-12_RKIND + + nCells = nCellsAll + nEdges = nEdgesAll + allocate(gradDensityEdge(nVertLevels, nEdges), & + dDensityDzTopOfCell(nVertLevels+1, nCells+1), & + gradDensityTopOfEdge(nVertLevels+1, nEdges), & + dDensityDzTopOfEdge(nVertLevels+1, nEdges), & + gradZMidEdge(nVertLevels, nEdges), & + gradZMidTopOfEdge(nVertLevels+1, nEdges)) + + !$omp parallel + !$omp do schedule(runtime) private(k) + do iEdge=1,nEdgesAll + do k=1,nVertLevels + gradDensityEdge(k,iEdge) = 0.0_RKIND + end do + end do + !$omp end do + !$omp end parallel + + nEdges = nEdgesHalo(2) + !$omp parallel + !$omp do schedule(runtime) private(k, rtmp) + do iCell = 1, nCells + 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 + ! 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:minLevelCell(iCell),iCell) = 0.0_RKIND + dDensityDzTopOfCell(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND + end do + !$omp end do + !$omp end parallel + + nEdges = nEdgesHalo( 2 ) + + ! Interpolate dDensityDzTopOfCell to edge and layer interface + !$omp parallel + !$omp do schedule(runtime) private(k, cell1, cell2) + do iEdge = 1, nEdges + 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)) + end do + end do + !$omp end do + !$omp end parallel + + ! Compute density gradient (gradDensityEdge) and gradient of zMid (gradZMidEdge) + ! along the constant coordinate surface. + ! The computed variables lives at edge and mid-layer depth + !$omp parallel + !$omp do schedule(runtime) private(cell1, cell2, k) + do iEdge = 1, nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,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 + end do + !$omp end do + + !$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. minLevelEdgeBot(iEdge)) then + do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) + h1 = layerThickEdgeMean(k-1,iEdge) + h2 = layerThickEdgeMean(k,iEdge) + ! Using second-order interpolation below + gradDensityTopOfEdge(k,iEdge) = (h2 * gradDensityEdge(k-1,iEdge) + h1 * & + gradDensityEdge(k,iEdge)) / (h1 + h2) + gradZMidTopOfEdge(k,iEdge) = (h2 * gradZMidEdge(k-1,iEdge) + h1 * & + gradZMidEdge(k,iEdge)) / (h1 + h2) + end do + + ! 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(minLevelEdgeBot(iEdge),iEdge) = gradDensityEdge(minLevelEdgeBot(iEdge),iEdge) + gradDensityTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradDensityEdge(maxLevelEdgeTop(iEdge),iEdge) + gradZMidTopOfEdge(minLevelEdgeBot(iEdge),iEdge) = gradZMidEdge(minLevelEdgeBot(iEdge),iEdge) + gradZMidTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradZMidEdge(maxLevelEdgeTop(iEdge),iEdge) + end if + end do + !$omp end do + + !$omp do schedule(runtime) private(k) + do iEdge = 1, nEdges + if (maxLevelEdgeTop(iEdge) .GE. minLevelEdgeBot(iEdge)) then + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1 + gradBuoyEddy(k,iEdge) = gravity*(gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & + * gradZMidTopOfEdge(k,iEdge))/rho_sw + end do + end if + end do + !$omp end do + !$omp end parallel + deallocate(gradDensityEdge, & + dDensityDzTopOfCell, & + gradDensityTopOfEdge, & + dDensityDzTopOfEdge, & + gradZMidEdge, & + gradZMidTopOfEdge) + + end subroutine ocn_eddy_compute_buoyancy_gradient + + !*********************************************************************** + ! + ! routine ocn_eddy_compute_mixed_layer_depth + ! + !> \brief Computes the density threshold based mixed layer depth + !> \details + !> Computes the density threshold based mixed layer depth. This + !> quantity is required by the submesoscale and redi parameterizations + !> The calculation is moved out of analysis members for this reason + !*********************************************************************** + + subroutine ocn_eddy_compute_mixed_layer_depth(statePool) + + type (mpas_pool_type), pointer, intent(in) :: statePool + + integer :: iEdge, nEdges, refIndex, cell1, cell2, nCells, k, iCell + real(kind=RKIND) :: dDenThres, den_ref_lev + real(kind=RKIND),dimension(:), allocatable :: depth + integer, dimension(:), pointer :: landIceMask + real (kind=RKIND), dimension(:,:), pointer :: layerThickness + logical :: found_den_mld + real (kind=RKIND) :: dV, dVp1, refDepth, coeffs(2), mldTemp + + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + nCells = nCellsAll + + allocate(depth(nVertLevels)) + refDepth = config_mld_reference_depth + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(depth, refIndex, coeffs, found_den_mld, k, den_ref_lev, dVp1, dV, mldTemp) + do iCell = 1,nCells + + found_den_mld = .false. + + depth(1) = 0.5_RKIND*layerThickness(1,iCell) + do k=2,maxLevelCell(iCell) + depth(k) = depth(k-1) + 0.5_RKIND*(layerThickness(k-1,iCell) + & + layerThickness(k ,iCell)) + end do + + do k=1, maxLevelCell(iCell)-1 + if(depth(k+1) >= refDepth) then + coeffs(2) = (potentialDensity(k+1,iCell) - potentialDensity(k,iCell)) & + / (depth(k+1) - depth(k)) + coeffs(1) = potentialDensity(k,iCell) - coeffs(2)*depth(k) + + den_ref_lev = coeffs(2)*refDepth + coeffs(1) + + refIndex=k + exit + end if + end do + + do k=refIndex,maxLevelCell(iCell)-1 + + if(.not. found_den_mld .and. (potentialDensity(k+1,iCell) - den_ref_lev) .ge. & + config_mixedLayerDepths_crit_dens_threshold) then + dVp1 = (potentialDensity(k+1,iCell) - den_ref_lev) + dV = (potentialDensity(k ,iCell) - den_ref_lev) + + coeffs(2) = (depth(k+1) - depth(k)) / (dVp1 - dV) + coeffs(1) = depth(k) - coeffs(2)*dV + + mldTemp = coeffs(2)*config_mixedLayerDepths_crit_dens_threshold + coeffs(1) + mldTemp=min(mldTemp,depth(k+1)) !make sure MLD isn't deeper than zMid(k+1) + dThreshMLD(iCell)=abs(max(mldTemp,depth(k))) !MLD should be deeper than zMid(k) + indMLD(iCell) = k+1 + found_den_mld = .true. + exit + end if + end do + + ! if the mixed layer depth is not found, it is set to the depth of the bottom most level + if(.not. found_den_mld) then + dThreshMLD(iCell) = depth(maxLevelCell(iCell)) + indMLD(iCell) = maxLevelCell(iCell) + end if + end do !iCell + !$omp end do + !$omp end parallel + + nEdges = nEdgesHalo(1) + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(cell1,cell2) + do iEdge=1,nEdgesAll + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + indMLDedge(iEdge) = min(indMLD(cell1),indMLD(cell2)) + enddo + !$omp end do + !$omp end parallel + deallocate(depth) + end subroutine ocn_eddy_compute_mixed_layer_depth + +end module ocn_eddy_parameterization_helpers diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index 0e6cdca08cec..6e968c908fda 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -16,6 +16,7 @@ module ocn_gm use ocn_constants use ocn_config use ocn_diagnostics_variables + use ocn_submesoscale_eddies implicit none private @@ -146,12 +147,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray, nEdgesArray - ! gradDensityEdge: Normal gradient of density - ! units: none - real(kind=RKIND), dimension(:, :), allocatable :: gradDensityEdge, gradDensityTopOfEdge, & - gradDensityConstZTopOfEdge, dDensityDzTopOfCell, dDensityDzTopOfEdge, & - gradZMidEdge, gradZMidTopOfEdge - type(mpas_pool_type), pointer :: tracersPool real(kind=RKIND), dimension(:, :, :), pointer :: activeTracers integer, pointer :: indexTemperature, indexSalinity @@ -196,19 +191,10 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & nCells = nCellsArray(size(nCellsArray)) nEdges = nEdgesArray(size(nEdgesArray)) - allocate(gradDensityEdge(nVertLevels, nEdges), & - dDensityDzTopOfCell(nVertLevels+1, nCells+1), & - gradDensityTopOfEdge(nVertLevels+1, nEdges), & - dDensityDzTopOfEdge(nVertLevels+1, nEdges), & - gradZMidEdge(nVertLevels, nEdges), & - gradZMidTopOfEdge(nVertLevels+1, nEdges), & - gradDensityConstZTopOfEdge(nVertLevels+1, nEdges)) - !$omp parallel !$omp do schedule(runtime) private(k) do iEdge = 1, nEdges do k = 1, nVertLevels - gradDensityEdge(k, iEdge) = 0.0_RKIND normalGMBolusVelocity(k, iEdge) = 0.0_RKIND end do end do @@ -364,6 +350,12 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & slopeTaperDown = 1.0_RKIND + slopeTaperFactor*(slopeTaperDown - 1.0_RKIND) sfcTaper = min(RediKappaSfcTaper(k, cell1), RediKappaSfcTaper(k, cell2)) + if(k < indMLDedge(iEdge)) then + sfcTaper = 0.0_RKIND + else + sfcTaper = 1.0_RKIND + end if + sfcTaperUp = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND) sfcTaperDown = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND) @@ -424,90 +416,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & !-------------------------------------------------------------------- if (config_use_GM) then - nEdges = nEdgesArray(3) - !$omp parallel - !$omp do schedule(runtime) private(k, rtmp) - do iCell = 1, nCells - 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 - - ! 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:minLevelCell(iCell),iCell) = 0.0_RKIND - dDensityDzTopOfCell(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND - end do - !$omp end do - !$omp end parallel - - nEdges = nEdgesArray( 3 ) - - ! Interpolate dDensityDzTopOfCell to edge and layer interface - !$omp parallel - !$omp do schedule(runtime) private(k, cell1, cell2) - do iEdge = 1, nEdges - 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)) - end do - end do - !$omp end do - !$omp end parallel - - ! Compute density gradient (gradDensityEdge) and gradient of zMid (gradZMidEdge) - ! along the constant coordinate surface. - ! The computed variables lives at edge and mid-layer depth - !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, k) - do iEdge = 1, nEdges - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,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 - end do - !$omp end do - - !$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. minLevelEdgeBot(iEdge)) then - do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge) - h1 = layerThickEdgeMean(k-1,iEdge) - h2 = layerThickEdgeMean(k,iEdge) - ! Using second-order interpolation below - gradDensityTopOfEdge(k,iEdge) = (h2 * gradDensityEdge(k-1,iEdge) + h1 * & - gradDensityEdge(k,iEdge)) / (h1 + h2) - gradZMidTopOfEdge(k,iEdge) = (h2 * gradZMidEdge(k-1,iEdge) + h1 * & - gradZMidEdge(k,iEdge)) / (h1 + h2) - end do - - ! 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(minLevelEdgeBot(iEdge),iEdge) = gradDensityEdge(minLevelEdgeBot(iEdge),iEdge) - gradDensityTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradDensityEdge(maxLevelEdgeTop(iEdge),iEdge) - gradZMidTopOfEdge(minLevelEdgeBot(iEdge),iEdge) = gradZMidEdge(minLevelEdgeBot(iEdge),iEdge) - gradZMidTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradZMidEdge(maxLevelEdgeTop(iEdge),iEdge) - end if - end do - !$omp end do - - !$omp do schedule(runtime) private(k) - do iEdge = 1, nEdges - 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) - end do - end if - end do - !$omp end do - !$omp end parallel nEdges = nEdgesArray(3) ! For config_GM_closure = 'N2_dependent' use a scaling to taper gmBolusKappa @@ -806,8 +714,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw & - *gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) ! Second to next to the last rows do k = minLevelEdgeBot(iEdge) + 2, maxLevelEdgeTop(iEdge) - 1 @@ -819,8 +726,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw & - *gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) end do ! Last row @@ -831,8 +737,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) tridiagB(k - 1) = -2.0_RKIND*cGMphaseSpeed(iEdge)**2/(layerThickEdgeMean(k - 1, iEdge) & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw & - *gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) ! Total number of rows N = maxLevelEdgeTop(iEdge) - minLevelEdgeBot(iEdge) @@ -865,13 +770,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & deallocate (tridiagC) end if !end config_use_GM - deallocate(gradDensityEdge, & - dDensityDzTopOfCell, & - gradDensityTopOfEdge, & - dDensityDzTopOfEdge, & - gradZMidEdge, & - gradZMidTopOfEdge, & - gradDensityConstZTopOfEdge) call mpas_timer_stop('gm bolus velocity') @@ -988,14 +886,7 @@ subroutine ocn_GM_init(domain, err)!{{{ end if if (config_Redi_use_surface_taper) then - if (config_AM_mixedLayerDepths_enable .and. config_AM_mixedLayerDepths_Dthreshold) then - sfcTaperFactor = 1.0_RKIND - else - call mpas_log_write('Redi Surface tapering requires MLD AM enabled with dThresh option selected.', & - MPAS_LOG_CRIT) - err = 1 - call mpas_dmpar_finalize(domain%dminfo) - end if + sfcTaperFactor = 1.0_RKIND else sfcTaperFactor = 0.0_RKIND end if diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F new file mode 100644 index 000000000000..8d0a1de912b5 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F @@ -0,0 +1,181 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! + +module ocn_submesoscale_eddies + + use mpas_pool_routines + use mpas_derived_types + use mpas_constants + use mpas_threading + use mpas_timer + + use ocn_constants + use ocn_config + use ocn_diagnostics_variables + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_submesoscale_compute_velocity, & + ocn_submesoscale_init + + real(kind=RKIND) :: tau, LfMin, Ce, dsMax + + + !*********************************************************************** + + contains + + !*********************************************************************** + ! + ! routine ocn_submesoscale_compute_velocity + ! + !> \brief Computes velocity from submesoscale eddies + !> \details + !> This routine implements the Fox-Kemper et al 2011 submesoscale + !> parameterization (https://doi.org/10.1016/j.ocemod.2010.09.002). + !> The transport velocity from the submesoscales is given by a + !> stream function which is given as + !> + !> Psi = C_e*Delta_S/L_f*H^2 int(bGrad_H)/sqrt(f^2+tau^-2)*mu(z) + !> + !> here C_e is a specified efficiency, Delta_s is taken as + !> min(dcEdge,dsMax), H is the mixed layer depth, L_f is the subgrid + !> frontal width that is scaled to the coarse grid, f is the + !> coriolis parameter, and tau is a time scale specified to prevent + !> the singularity near the equator. the horizontal buoyancy gradient + !> is integrated over the mixed layer depth and mu(z) is a non-dimensional + !> shape function that is zero below the mixed layer depth + !> + !*********************************************************************** + + subroutine ocn_submesoscale_compute_velocity() + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2 + + real(kind=RKIND),dimension(:),allocatable :: & + streamFunction, & + zEdge + + real(kind=RKIND) :: mu, zMLD, bvfML, hML, gradBuoyML, ds, Lf, bvfAv, & + hAv, bldEdge, ustarEdge + + nEdges = nEdgesHalo(2) + + allocate(streamFunction(nVertLevels+1)) + allocate(zEdge(nVertLevels+1)) + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(cell1,cell2,bvfML,zEdge,streamFunction,gradBuoyML,hML,bvfAv) + !$omp private(bvfML,zMLD,Lf,ds,mu,k) + do iEdge=1,nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + zEdge(:) = 0.0_RKIND + streamFunction(:) = 0.0_RKIND + bvfML = 0.0_RKIND + gradBuoyML = 0.0_RKIND + bvfAv = sqrt(0.5_RKIND*(max(BruntVaisalaFreqTop(minLevelCell(cell1)+1,cell1),1.0E-20_RKIND) + & + max(BruntVaisalaFreqTop(minLevelCell(cell2)+1,cell2),1.0E-20_RKIND))) + hML = 0.5_RKIND*layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge) + bvfML = bvfML + hML*bvfAv + if (minLevelEdgeTop(iEdge) .ge. 1) then + gradBuoyML = hML*gradBuoyEddy(minLevelEdgeTop(iEdge)+1,iEdge) + bvfML = hML*bvfAv + else + cycle + end if + + do k = minLevelEdgeBot(iEdge)+2,indMLDedge(iEdge) + hAv = 0.5_RKIND*(layerThickEdgeMean(k,iEdge) + layerThickEdgeMean(k-1,iEdge)) + bvfML = bvfML + hAv*bvfAv + hML = hML + hAv + bvfAv = sqrt(0.5_RKIND*(max(BruntVaisalaFreqTop(k,cell1),1.0E-20_RKIND) + & + max(BruntVaisalaFreqTop(k,cell2),1.0E-20_RKIND))) + gradBuoyML = gradBuoyML + hAv*gradBuoyEddy(k,iEdge) + end do + bvfML = bvfML / (1.0E-20_RKIND + hML) + gradBuoyML = gradBuoyML / (1.0E-20_RKIND + hML) + + !compute depths and shape function + + do k = minLevelEdgeTop(iEdge)+1,maxLevelEdgeTop(iEdge)+1 + zEdge(k) = zedge(k-1) - layerThickEdgeMean(k-1,iEdge) + end do + + zMLD = 0.5_RKIND*(dThreshMLD(cell1)+dThreshMLD(cell2)) + Lf = max(Lfmin, abs(gradBuoyML)*zMLD / (1.0E-20_RKIND + fEdge(iEdge)**2), & + bvfML*zMLD / abs(fEdge(iEdge))) + + ds = min(dcEdge(iEdge),dsMax) + + do k = minLevelEdgeTop(iEdge)+1,maxLevelEdgeTop(iEdge) + mu = max(0.0_RKIND,(1.0_RKIND - (2.0_RKIND*zEdge(k) / zMLD + 1.0_RKIND)**2.0)* & + (1.0_RKIND + 5.0_RKIND/21.0_RKIND*(2.0_RKIND*zEdge(k)/zMLD + 1.0_RKIND)**2.0)) + streamFunction(k) = Ce*ds/Lf*zMLD**2.0*gradBuoyML/sqrt(fEdge(iEdge)**2.0 + tau**(-2.0))*mu + end do + + ! integrate in vertical to get the velocity + do k = minLevelEdgeTop(iEdge),maxLevelEdgeTop(iEdge) + normalMLEvelocity(k,iEdge) = -(streamFunction(k) - streamFunction(k+1)) / layerThickEdgeMean(k,iEdge) + end do + + end do!iEdge loop + !$omp end do + !$omp end parallel + + deallocate(streamFunction) + deallocate(zEdge) + + end subroutine ocn_submesoscale_compute_velocity + + !*********************************************************************** + ! + ! routine ocn_submesoscale_init + ! + !> \brief Submesoscale parameterization init + !> \details + !> Initializes parameters related to the submesoscale eddy parameterization + !> + !*********************************************************************** + + subroutine ocn_submesoscale_init(err) + + integer, intent(out) :: err !< Output: error flag + + err = 0 + Ce = config_submesoscale_ce + Lfmin = config_submesoscale_lfmin + dsmax = config_submesoscale_ds_max + tau = config_submesoscale_tau + + end subroutine ocn_submesoscale_init + +end module ocn_submesoscale_eddies diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F index 47d17483e276..b6a558b4a522 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -238,14 +238,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea 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)) + RediKappaScaling(k + 1, cell1) + RediKappaScaling(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) + & - RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2)) + RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))!* & do iTr = 1, nTracers tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) @@ -280,9 +276,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea 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* & - (RediKappaSfcTaper(k, cell1) + RediKappaSfcTaper(k, cell2) + & - RediKappaSfcTaper(k + 1, cell1) + RediKappaSfcTaper(k + 1, cell2)) + RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))!*0.25_RKIND* & do iTr = 1, nTracers ! \kappa_2 \nabla \phi on edge @@ -319,9 +313,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea k = maxLevelEdgeTop(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)) + RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))!* & do iTr = 1, nTracers tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) @@ -376,8 +368,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea ! 2.0 in next line is because a dot product on a C-grid ! requires a factor of 1/2 to average to the cell center. flux_term3 = rediKappaCell(iCell)*2.0_RKIND* & - (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) & - * invAreaCell - RediKappaSfcTaper(k + 1, iCell)*RediKappaScaling(k + 1, iCell) & + (RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) & + * invAreaCell - RediKappaScaling(k + 1, iCell) & *fluxRediZTop(iTr, k + 1) * invAreaCell) redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3 @@ -446,9 +438,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* & - (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)* & - redi_term3_topOfCell(iTr, k, iCell) - RediKappaSfcTaper(k + 1, iCell) & - *RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell)) + (RediKappaScaling(k, iCell)* & + redi_term3_topOfCell(iTr, k, iCell) - & + RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell)) end do end do From d39413ded9fceda62552a80441d0788681e003fa Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Wed, 27 Jul 2022 15:57:02 -0500 Subject: [PATCH 02/28] Fixes CF compliance on new variables --- cime | 2 +- components/elm/src/external_models/fates | 2 +- components/elm/src/external_models/sbetr | 2 +- components/mpas-ocean/src/Registry.xml | 22 +++++++++++----------- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/cime b/cime index a00a6ed2fb2e..640f471c02f7 160000 --- a/cime +++ b/cime @@ -1 +1 @@ -Subproject commit a00a6ed2fb2e8a8dceedc3973585d7cc3d867707 +Subproject commit 640f471c02f7663497c652aedf39a6564e7e800d diff --git a/components/elm/src/external_models/fates b/components/elm/src/external_models/fates index c11c448bcb25..fccbd303069b 160000 --- a/components/elm/src/external_models/fates +++ b/components/elm/src/external_models/fates @@ -1 +1 @@ -Subproject commit c11c448bcb2531220134cd3da70cac6712a5b759 +Subproject commit fccbd303069bf66895903ce7e145fce6e1d54af6 diff --git a/components/elm/src/external_models/sbetr b/components/elm/src/external_models/sbetr index 13fff9208624..1934b9ec58d2 160000 --- a/components/elm/src/external_models/sbetr +++ b/components/elm/src/external_models/sbetr @@ -1 +1 @@ -Subproject commit 13fff9208624ca5e7da1094f0d93043f8cd58926 +Subproject commit 1934b9ec58d2b82db4d52e086e5e5fe31a632d23 diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 1515ad1b2f59..076c7dde2754 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -344,7 +344,7 @@ /> - @@ -352,7 +352,7 @@ description="timescale for frictional slumping of front (in seconds)" possible_values="positive reals, between 1-10 days" /> - @@ -2703,7 +2703,7 @@ description="vertical tracer-transport velocity defined at center (horizontally) and top (vertically) of cell. This is not the vertical ALE transport, but is Eulerian (fixed-frame) in the vertical, and computed from the continuity equation from the horizontal GM Bolus velocity." packages="forwardMode;analysisMode" /> - - @@ -2865,27 +2865,27 @@ description="count of times redi limiter is invoked on a timestep" packages="gm" /> - - - - - - @@ -2986,7 +2986,7 @@ - Date: Wed, 27 Jul 2022 16:12:05 -0500 Subject: [PATCH 03/28] Fixes space vs tabs in Registry --- components/mpas-ocean/src/Registry.xml | 34 +++++++++++++------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 076c7dde2754..e00dbec531cf 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -344,24 +344,24 @@ /> - - - - - + @@ -2838,8 +2838,8 @@ description="Bolus velocity in Gent-McWilliams eddy parameterization" packages="forwardMode;analysisMode" /> - - - + Date: Wed, 27 Jul 2022 22:17:30 -0500 Subject: [PATCH 04/28] Fixes cime commit --- cime | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cime b/cime index 640f471c02f7..a00a6ed2fb2e 160000 --- a/cime +++ b/cime @@ -1 +1 @@ -Subproject commit 640f471c02f7663497c652aedf39a6564e7e800d +Subproject commit a00a6ed2fb2e8a8dceedc3973585d7cc3d867707 From 04d0086f753dfc4d35b7bd2f451c7db151ad8365 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Wed, 27 Jul 2022 22:21:07 -0500 Subject: [PATCH 05/28] Fixes ELM commits --- components/elm/src/external_models/fates | 2 +- components/elm/src/external_models/sbetr | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/components/elm/src/external_models/fates b/components/elm/src/external_models/fates index fccbd303069b..c11c448bcb25 160000 --- a/components/elm/src/external_models/fates +++ b/components/elm/src/external_models/fates @@ -1 +1 @@ -Subproject commit fccbd303069bf66895903ce7e145fce6e1d54af6 +Subproject commit c11c448bcb2531220134cd3da70cac6712a5b759 diff --git a/components/elm/src/external_models/sbetr b/components/elm/src/external_models/sbetr index 1934b9ec58d2..13fff9208624 160000 --- a/components/elm/src/external_models/sbetr +++ b/components/elm/src/external_models/sbetr @@ -1 +1 @@ -Subproject commit 1934b9ec58d2b82db4d52e086e5e5fe31a632d23 +Subproject commit 13fff9208624ca5e7da1094f0d93043f8cd58926 From 269e8d206721d592a0dd8126f9bebcf9a9d0d08c Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 4 Aug 2022 09:55:45 -0500 Subject: [PATCH 06/28] Fix to OMP directive --- .../mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F index 8d0a1de912b5..3409398e4ff2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F @@ -92,8 +92,8 @@ subroutine ocn_submesoscale_compute_velocity() !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1,cell2,bvfML,zEdge,streamFunction,gradBuoyML,hML,bvfAv) - !$omp private(bvfML,zMLD,Lf,ds,mu,k) + !$omp private(cell1,cell2,bvfML,zEdge,streamFunction,gradBuoyML,hML,bvfAv, & + !$omp zMLD,Lf,ds,mu,k) do iEdge=1,nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) From ff94d3b2863b6a099da2fc54fb82e854b8a41a6a Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Sun, 7 Aug 2022 22:42:23 -0500 Subject: [PATCH 07/28] Reverts MOC calculation and disables mle param --- .../bld/namelist_files/namelist_defaults_mpaso.xml | 2 +- .../src/analysis_members/mpas_ocn_moc_streamfunction.F | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 35679ddd12f1..6df175e82fda 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -153,7 +153,7 @@ 40e3 -.true. +.false. 172800 0.06 1000.0 diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F index 5a2f71d71995..3daa9c6cdd91 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F @@ -507,12 +507,10 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ * 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2)) end do end do - mocStreamValLatAndDepthRegionLocal(1, nVertLevels, iTransect) = & - - sumTransport(nVertLevels, iTransect) - do k = nVertLevels-1,1,-1 + do k = 2, nVertLevels mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & - - sumTransport(k + 1, iTransect) + + sumTransport(k + 1, iTransect) end do end do From bcaf4bca1b674025cf7336a000c159924970b6e6 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Sun, 7 Aug 2022 23:07:36 -0500 Subject: [PATCH 08/28] Adds submesoscale induced MOC calculation --- components/mpas-ocean/cime_config/buildnml | 2 + .../Registry_moc_streamfunction.xml | 8 ++ .../mpas_ocn_moc_streamfunction.F | 122 ++++++++++++++++-- 3 files changed, 121 insertions(+), 11 deletions(-) diff --git a/components/mpas-ocean/cime_config/buildnml b/components/mpas-ocean/cime_config/buildnml index 2e04d3516cf4..343ce7c6615a 100755 --- a/components/mpas-ocean/cime_config/buildnml +++ b/components/mpas-ocean/cime_config/buildnml @@ -1125,6 +1125,8 @@ def buildnml(case, caseroot, compname): lines.append(' ') lines.append(' ') + lines.append(' ') + lines.append(' ') lines.append(' ') lines.append(' ') lines.append(' ') diff --git a/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml b/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml index 9019ebce28c2..a2b32ebe6ed0 100644 --- a/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml +++ b/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml @@ -61,9 +61,15 @@ + + @@ -89,6 +95,8 @@ + + diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F index 3daa9c6cdd91..5a5f3f538329 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F @@ -578,7 +578,6 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ call mpas_pool_get_dimension(block % dimensions, 'nEdgesSolve', nEdgesSolve) call mpas_pool_get_array(transectPool,'transectEdgeMaskSigns',transectEdgeMaskSigns) call mpas_pool_get_array(transectPool,'transectEdgeMasks',transectEdgeMasks) - normalVelocity => normalGMBolusVelocity call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) @@ -602,18 +601,14 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ sumTransport(k,iTransect) = sumTransport(k,iTransect) + & transectEdgeMaskSigns(currentTransect,iEdge) & * transectEdgeMasks(currentTransect, iEdge) & - * (normalMLEvelocity(k,iEdge) + normalVelocity(k,iEdge)) & - *dvEdge(iEdge) & + * normalGMBolusVelocity(k,iEdge)*dvEdge(iEdge) & * 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2)) end do end do - ! do k = 2, nVertLevels - - mocStreamValLatAndDepthRegionLocal(1, nVertLevels, iTransect) = -sumTransport(nVertLevels,iTransect) - do k = nVertLevels-1,1,-1 + do k = 2, nVertLevels mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & - - sumTransport(k + 1, iTransect) + + sumTransport(k + 1, iTransect) end do end do @@ -625,12 +620,11 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ do i = 1, regionsInAddGroup currentRegion = regionsInGroup(i, regionGroupNumber) sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + & - ((vertGMBolusVelocityTop(k, iCell) +vertMLEBolusVelocityTop(k,iCell)) * & + (vertGMBolusVelocityTop(k, iCell) * & areaCell(iCell) * regionCellMasks(currentRegion, iCell)) end do sumVertBinVelocity(iBin, k) = sumVertBinVelocity(iBin, k) + & - ((vertGMBolusVelocityTop(k, iCell) + vertMLEBolusVelocityTop(k,iCell)) & - * areaCell(iCell)) + (vertGMBolusVelocityTop(k, iCell) * areaCell(iCell)) end do end do @@ -667,6 +661,112 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ endif !config_use_GM + if(config_submesoscale_enable) THEN !compute submesoscale eddy bolus contribution to the streamfunction + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'mocStreamfunctionAM', mocStreamfunctionAMPool) + + call mpas_pool_get_array(mocStreamfunctionAMPool, 'binBoundaryMocStreamfunction', binBoundaryMocStreamfunction) + + binWidth = (binBoundaryMocStreamfunction(nMocStreamfunctionBinsUsed + 1) - binBoundaryMocStreamfunction(1)) & + / nMocStreamfunctionBinsUsed + + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) + + call mpas_pool_get_array(meshPool, 'latCell', latCell) + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + + call mpas_pool_get_array(regionPool, 'regionCellMasks', regionCellMasks) + + !!!! TRANSECT DOMAINSPLIT VARIABLES + call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) + call mpas_pool_get_dimension(block % dimensions, 'nEdgesSolve', nEdgesSolve) + call mpas_pool_get_array(transectPool,'transectEdgeMaskSigns',transectEdgeMaskSigns) + call mpas_pool_get_array(transectPool,'transectEdgeMasks',transectEdgeMasks) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) + + !!!! TRANSECT CALCULATION + sumTransport = 0.0_RKIND + mocStreamValLatAndDepthLocal = 0.0_RKIND + mocStreamValLatAndDepthTotal = 0.0_RKIND + sumVertBinVelocity = 0.0_RKIND + mocStreamValLatAndDepthRegionLocal = 0.0_RKIND + mocStreamValLatAndDepthRegionTotal = 0.0_RKIND + sumVertBinVelocityRegion = 0.0_RKIND + + do iTransect = 1,transectsInAddGroup + currentTransect = transectsInGroup(iTransect, transectGroupNumber) + do iEdge = 1,nEdgesSolve + c1 = cellsOnEdge(1,iEdge) + c2 = cellsOnEdge(2,iEdge) + do k = 1, maxLevelEdgeTop(iEdge) + sumTransport(k,iTransect) = sumTransport(k,iTransect) + & + transectEdgeMaskSigns(currentTransect,iEdge) & + * transectEdgeMasks(currentTransect, iEdge) & + * normalMLEvelocity(k,iEdge)*dvEdge(iEdge) & + * 0.5_RKIND*(layerThickness(k,c1) + layerThickness(k,c2)) + end do + end do + do k = 2, nVertLevels + mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & + mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & + + sumTransport(k + 1, iTransect) + end do + end do + + !!!! END TRANSECT CALCULATION + + do iCell = 1,nCellsSolve + iBin = MAX(int((latCell(iCell) - binBoundaryMocStreamfunction(1)) / binWidth), 2) + do k = 1,maxLevelCell(iCell) + do i = 1, regionsInAddGroup + currentRegion = regionsInGroup(i, regionGroupNumber) + sumVertBinVelocityRegion(iBin, k, i) = sumVertBinVelocityRegion(iBin, k, i) + & + (vertMLEBolusVelocityTop(k,iCell) * & + areaCell(iCell) * regionCellMasks(currentRegion, iCell)) + end do + sumVertBinVelocity(iBin, k) = sumVertBinVelocity(iBin, k) + & + (vertMLEBolusVelocityTop(k,iCell) * areaCell(iCell)) + end do + end do + + do i = 1, regionsInAddGroup + do k = 1,nVertLevels + do iBin = 2, nMocStreamfunctionBinsUsed + 1 + mocStreamValLatAndDepthLocal(iBin, k) = mocStreamValLatAndDepthLocal(iBin-1, k) & + + sumVertBinVelocity(iBin, k) + mocStreamValLatAndDepthRegionLocal(iBin, k, i) = mocStreamValLatAndDepthRegionLocal(iBin-1, k, i) & + + sumVertBinVelocityRegion(iBin, k, i) + end do + end do + end do + + block => block % next + end do + + call mpas_dmpar_sum_real_array(dminfo, nVertLevels * (nMocStreamfunctionBinsUsed + 1), mocStreamValLatAndDepthLocal, & + mocStreamvalLatAndDepthTotal) + + call mpas_dmpar_sum_real_array(dminfo, nVertLevels * (nMocStreamfunctionBinsUsed + 1) * maxRegionsInGroup, & + mocStreamValLatAndDepthRegionLocal, mocStreamvalLatAndDepthRegionTotal) + + call mpas_pool_get_subpool(domain % blocklist % structs, 'mocStreamfunctionAM', mocStreamfunctionAMPool) + call mpas_pool_get_array(mocStreamfunctionAMPool, 'mocStreamvalLatAndDepthMLE', mocStreamvalLatAndDepthMLE) + mocStreamvalLatAndDepthMLE = mocStreamvalLatAndDepthTotal * m3ps_to_Sv + + call mpas_pool_get_array(mocStreamfunctionAMPool, 'mocStreamvalLatAndDepthRegionMLE', mocStreamvalLatAndDepthRegionMLE) + mocStreamvalLatAndDepthRegionMLE = mocStreamvalLatAndDepthRegionTotal * m3ps_to_Sv + + !Add submesoscale eddy bolus contribution to resolved streamfunction to create total streamfunction + mocStreamvalLatAndDepthRegion = mocStreamvalLatAndDepthRegion + mocStreamvalLatAndDepthRegionMLE + mocStreamvalLatAndDepth = mocStreamvalLatAndDepth + mocStreamvalLatAndDepthMLE + + endif !config_use_GM + deallocate(mocStreamvalLatAndDepthTotal) deallocate(mocStreamvalLatAndDepthLocal) deallocate(sumVertBinVelocity) From fadff8fa1aa40f5fe5baf060f07e2ff6e034569a Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Sun, 7 Aug 2022 23:28:58 -0500 Subject: [PATCH 09/28] Removes some extra variables --- components/mpas-ocean/src/Registry.xml | 27 ++------- .../mpas_ocn_time_integration_rk4.F | 4 +- .../mpas_ocn_time_integration_si.F | 4 +- .../mpas_ocn_time_integration_split.F | 4 +- .../src/shared/mpas_ocn_diagnostics.F | 24 ++++---- .../shared/mpas_ocn_diagnostics_variables.F | 57 +++++++------------ .../mpas_ocn_eddy_parameterization_helpers.F | 13 +---- .../src/shared/mpas_ocn_submesoscale_eddies.F | 5 +- 8 files changed, 47 insertions(+), 91 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index e00dbec531cf..36368086850e 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2869,18 +2869,6 @@ description="horizontal gradient of buoyancy at cell edge and interfaces in vertical. this is used for the GM parameterization and the submesoscale eddy parameterization" packages="forwardMode" /> - - - - - - - diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 454ae89ac1eb..2bd33ab53476 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -819,8 +819,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call ocn_time_average_coupled_accumulate(statePool, forcingPool, 2) - if (config_use_GM) then - call ocn_reconstruct_gm_vectors(meshPool) + if (config_use_GM .or. config_submesoscale_enable) then + call ocn_reconstruct_eddy_vectors(meshPool) end if block => block % next diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index cda0954fd1ef..f8b39b21ab96 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -2566,8 +2566,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ #endif call ocn_time_average_coupled_accumulate(statePool,forcingPool,2) - if (config_use_GM) then - call ocn_reconstruct_gm_vectors(meshPool) + if (config_use_GM .or. config_submesoscale_enable) then + call ocn_reconstruct_eddy_vectors(meshPool) end if if (trim(config_land_ice_flux_mode) == 'coupled') then diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 556a73dbd6b2..7a19a5db9530 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -2630,8 +2630,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call ocn_time_average_coupled_accumulate(statePool,forcingPool,2) - if (config_use_GM) then - call ocn_reconstruct_gm_vectors(meshPool) + if (config_use_GM .or. config_submesoscale_enable) then + call ocn_reconstruct_eddy_vectors(meshPool) end if if (trim(config_land_ice_flux_mode) == 'coupled') then diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 1e79bd9d6cdc..cb0c3bb77e79 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -61,7 +61,7 @@ module ocn_diagnostics ocn_vert_transport_velocity_top, & ocn_fuperp, & ocn_filter_btr_mode_tend_vel, & - ocn_reconstruct_gm_vectors, & + ocn_reconstruct_eddy_vectors, & ocn_compute_kpp_input_fields, & ocn_validate_state, & ocn_build_log_filename, & @@ -3095,7 +3095,7 @@ end subroutine ocn_compute_land_ice_flux_input_fields!}}} !*********************************************************************** ! -! routine ocn_reconstruct_gm_vectors +! routine ocn_reconstruct_eddy_vectors ! !> \brief Computes cell-centered vector diagnostics !> \author Mark Petersen @@ -3105,11 +3105,11 @@ end subroutine ocn_compute_land_ice_flux_input_fields!}}} ! !----------------------------------------------------------------------- - subroutine ocn_reconstruct_gm_vectors(meshPool) !{{{ + subroutine ocn_reconstruct_eddy_vectors(meshPool) !{{{ type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information - call mpas_timer_start('reconstruct gm vecs') + call mpas_timer_start('reconstruct eddy vecs') call mpas_reconstruct(meshPool, normalTransportVelocity, & transportVelocityX, & @@ -3120,17 +3120,17 @@ subroutine ocn_reconstruct_gm_vectors(meshPool) !{{{ ) call mpas_reconstruct(meshPool, normalGMBolusVelocity, & - GMBolusVelocityX, & - GMBolusVelocityY, & - GMBolusVelocityZ, & + eddyVelocityX, & + eddyVelocityY, & + eddyVelocityZ, & GMBolusVelocityZonal, & GMBolusVelocityMeridional & ) call mpas_reconstruct(meshPool, normalMLEvelocity, & - mleVelocityX, & - mleVelocityY, & - mleVelocityZ, & + eddyVelocityX, & + eddyVelocityY, & + eddyVelocityZ, & mleVelocityZonal, & mleVelocityMeridional & ) @@ -3143,9 +3143,9 @@ subroutine ocn_reconstruct_gm_vectors(meshPool) !{{{ GMStreamFuncMeridional & ) - call mpas_timer_stop('reconstruct gm vecs') + call mpas_timer_stop('reconstruct eddy vecs') - end subroutine ocn_reconstruct_gm_vectors!}}} + end subroutine ocn_reconstruct_eddy_vectors!}}} !*********************************************************************** diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index 644257970ff5..ad84e5c92308 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -99,10 +99,10 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:), pointer :: & transportVelocityX, transportVelocityY, transportVelocityZ, transportVelocityZonal, & - transportVelocityMeridional, GMBolusVelocityX, GMBolusVelocityY, & - GMBolusVelocityZ, GMBolusVelocityZonal, GMBolusVelocityMeridional, gmStreamFuncTopOfEdge, & + transportVelocityMeridional, eddyVelocityX, eddyVelocityY, & + eddyVelocityZ, GMBolusVelocityZonal, GMBolusVelocityMeridional, gmStreamFuncTopOfEdge, & GMStreamFuncX, GMStreamFuncY, GMStreamFuncZ, GMStreamFuncZonal, GMStreamFuncMeridional, & - mleVelocityX, mleVelocityY, mleVelocityZ, mleVelocityZonal, mleVelocityMeridional + mleVelocityZonal, mleVelocityMeridional real (kind=RKIND), dimension(:,:), pointer :: rx1Edge real (kind=RKIND), dimension(:,:), pointer :: rx1Cell @@ -127,7 +127,7 @@ module ocn_diagnostics_variables real(kind=RKIND), dimension(:), pointer :: cGMphaseSpeed real(kind=RKIND), dimension(:,:,:), pointer :: slopeTriadUp, slopeTriadDown - integer, dimension(:), pointer :: indMLD, indMLDedge + integer, dimension(:), pointer :: indMLD real(kind=RKIND), dimension(:), pointer :: dThreshMLD real (kind=RKIND), dimension(:,:), pointer :: velocityX, velocityY, velocityZ @@ -417,12 +417,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e GMStreamFuncZonal) call mpas_pool_get_array(diagnosticsPool, 'GMStreamFuncMeridional', & GMStreamFuncMeridional) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityX', & - mleVelocityX) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityY', & - mleVelocityY) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZ', & - mleVelocityZ) call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZonal', & mleVelocityZonal) call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & @@ -437,12 +431,12 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e mleVelocityZonal) call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & mleVelocityMeridional) - call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityX', & - GMBolusVelocityX) - call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityY', & - GMBolusVelocityY) - call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityZ', & - GMBolusVelocityZ) + call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityX', & + eddyVelocityX) + call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityY', & + eddyVelocityY) + call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityZ', & + eddyVelocityZ) call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityZonal', & GMBolusVelocityZonal) call mpas_pool_get_array(diagnosticsPool, 'GMBolusVelocityMeridional', & @@ -623,7 +617,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e 'modifyLandIcePressureMask', & modifyLandIcePressureMask) call mpas_pool_get_array(diagnosticsPool, 'indMLD', indMLD) - call mpas_pool_get_array(diagnosticsPool, 'indMLDedge', indMLDedge) call mpas_pool_get_array(diagnosticsPool, 'dThreshMLD', dThreshMLD) call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', & surfacePressure) @@ -724,12 +717,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & - !$acc mleVelocityX, & - !$acc mleVelocityY, & - !$acc mleVelocityZ, & - !$acc GMBolusVelocityX, & - !$acc GMBolusVelocityY, & - !$acc GMBolusVelocityZ, & + !$acc eddyVelocityX, & + !$acc eddyVelocityY, & + !$acc eddyVelocityZ, & !$acc vertNonLocalFlux, & !$acc vertViscTopOfEdge, & !$acc gradSSHMeridional, & @@ -859,7 +849,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc zTop, & !$acc xtime, & !$acc indMLD, & - !$acc indMLDedge, & !$acc dThreshMLD, & !$acc density, & !$acc betaEdge, & @@ -976,12 +965,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & - !$acc mleVelocityX, & - !$acc mleVelocityY, & - !$acc mleVelocityZ, & - !$acc GMBolusVelocityX, & - !$acc GMBolusVelocityY, & - !$acc GMBolusVelocityZ, & + !$acc eddyVelocityX, & + !$acc eddyVelocityY, & + !$acc eddyVelocityZ, & !$acc vertNonLocalFlux, & !$acc vertViscTopOfEdge, & !$acc gradSSHMeridional, & @@ -1111,7 +1097,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc zTop, & !$acc xtime, & !$acc indMLD, & - !$acc indMLDedge, & !$acc dThreshMLD, & !$acc density, & !$acc betaEdge, & @@ -1186,12 +1171,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ normRelVortEdge, & surfaceVelocity, & vertVelocityTop, & - GMBolusVelocityX, & - GMBolusVelocityY, & - GMBolusVelocityZ, & - mleVelocityX, & - mleVelocityY, & - mleVelocityZ, & + eddyVelocityX, & + eddyVelocityY, & + eddyVelocityZ, & vertNonLocalFlux, & vertViscTopOfEdge, & gradSSHMeridional, & @@ -1312,7 +1294,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ zTop, & xtime, & indMLD, & - indMLDedge, & dThreshMLD, & density, & betaEdge, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F index 5c1f12faedba..d8edd70eaa5e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F @@ -62,7 +62,7 @@ subroutine ocn_eddy_compute_buoyancy_gradient() ! !----------------------------------------------------------------- - integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2, indMLDedge + integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2 ! gradDensityEdge: Normal gradient of density ! units: none real(kind=RKIND), dimension(:, :), allocatable :: gradDensityEdge, gradDensityTopOfEdge, & @@ -265,17 +265,6 @@ subroutine ocn_eddy_compute_mixed_layer_depth(statePool) !$omp end do !$omp end parallel - nEdges = nEdgesHalo(1) - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(cell1,cell2) - do iEdge=1,nEdgesAll - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - indMLDedge(iEdge) = min(indMLD(cell1),indMLD(cell2)) - enddo - !$omp end do - !$omp end parallel deallocate(depth) end subroutine ocn_eddy_compute_mixed_layer_depth diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F index 3409398e4ff2..1dbf51b94fd4 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F @@ -76,7 +76,7 @@ subroutine ocn_submesoscale_compute_velocity() ! !----------------------------------------------------------------- - integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2 + integer :: k, nEdges, nCells, iCell, iEdge, cell1, cell2, indMLDedge real(kind=RKIND),dimension(:),allocatable :: & streamFunction, & @@ -113,7 +113,8 @@ subroutine ocn_submesoscale_compute_velocity() cycle end if - do k = minLevelEdgeBot(iEdge)+2,indMLDedge(iEdge) + indMLDedge = min(indMLD(cell1),indMLD(cell2)) + do k = minLevelEdgeBot(iEdge)+2,indMLDedge hAv = 0.5_RKIND*(layerThickEdgeMean(k,iEdge) + layerThickEdgeMean(k-1,iEdge)) bvfML = bvfML + hAv*bvfAv hML = hML + hAv From d26f43b9f453138ae07728c68ec92eb8890dc849 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Sun, 7 Aug 2022 23:54:15 -0500 Subject: [PATCH 10/28] Redi cleanup --- components/mpas-ocean/src/Registry.xml | 17 +++++++++-------- .../src/driver/mpas_ocn_core_interface.F | 11 +++++++++++ .../src/shared/mpas_ocn_tracer_hmix_redi.F | 4 ++-- 3 files changed, 22 insertions(+), 10 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 36368086850e..d631aeaf1010 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1519,6 +1519,7 @@ + @@ -2867,35 +2868,35 @@ /> Date: Mon, 8 Aug 2022 09:04:43 -0500 Subject: [PATCH 11/28] Reverts PE compset changes for chrysalis --- .../mpas-ocean/cime_config/config_pes.xml | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/components/mpas-ocean/cime_config/config_pes.xml b/components/mpas-ocean/cime_config/config_pes.xml index dcdf96722438..7cd87daed67d 100644 --- a/components/mpas-ocean/cime_config/config_pes.xml +++ b/components/mpas-ocean/cime_config/config_pes.xml @@ -235,34 +235,34 @@ mpas-ocean+chrysalis: standard-res, compset=DATM+MPASO, 15x32x2 NODESxMPIxOMP - 64 + 32 64 - 320 - 320 - 320 - 320 - 640 - 320 - 320 - 320 + 160 + 160 + 160 + 160 + 320 + 160 + 160 + 160 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 0 0 0 0 - 320 + 160 0 0 0 From 12c930660368a8c6f4d9d66e2940c00e5857bbf2 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 8 Aug 2022 23:20:04 -0500 Subject: [PATCH 12/28] Moves vertGM, vertMLE to packages --- components/mpas-ocean/driver/ocn_comp_mct.F | 2 +- components/mpas-ocean/src/Registry.xml | 11 +- .../src/shared/mpas_ocn_diagnostics.F | 111 +++++++++++++++--- 3 files changed, 99 insertions(+), 25 deletions(-) diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index 3517de9bafbb..2b9eb424a7e2 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -974,7 +974,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr) call ocn_eddy_compute_mixed_layer_depth(statePool) - if (config_use_GM.or.config_submesoscale_enable) then + if (config_use_GM .or. config_submesoscale_enable) then call ocn_eddy_compute_buoyancy_gradient() end if diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index d631aeaf1010..d3143ba56fa7 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2702,12 +2702,11 @@ /> Date: Mon, 8 Aug 2022 23:23:11 -0500 Subject: [PATCH 13/28] Fixes alignment --- .../shared/mpas_ocn_diagnostics_variables.F | 42 +++++++------------ 1 file changed, 16 insertions(+), 26 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index ad84e5c92308..3850155fb4a3 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -421,17 +421,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e mleVelocityZonal) call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & mleVelocityMeridional) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityX', & - mleVelocityX) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityY', & - mleVelocityY) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZ', & - mleVelocityZ) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityZonal', & - mleVelocityZonal) - call mpas_pool_get_array(diagnosticsPool, 'mleVelocityMeridional', & - mleVelocityMeridional) - call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityX', & + call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityX', & eddyVelocityX) call mpas_pool_get_array(diagnosticsPool, 'eddyVelocityY', & eddyVelocityY) @@ -717,9 +707,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & - !$acc eddyVelocityX, & - !$acc eddyVelocityY, & - !$acc eddyVelocityZ, & + !$acc eddyVelocityX, & + !$acc eddyVelocityY, & + !$acc eddyVelocityZ, & !$acc vertNonLocalFlux, & !$acc vertViscTopOfEdge, & !$acc gradSSHMeridional, & @@ -739,7 +729,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc BruntVaisalaFreqTop, & !$acc vertAleTransportTop, & !$acc GMBolusVelocityZonal, & - !$acc mleVelocityZonal, & + !$acc mleVelocityZonal, & !$acc bulkRichardsonNumber, & !$acc gmStreamFuncTopOfEdge, & !$acc indexSSHGradientZonal, & @@ -760,7 +750,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc bulkRichardsonNumberBuoy, & !$acc tracersSurfaceLayerValue, & !$acc GMBolusVelocityMeridional, & - !$acc mleVelocityMeridional, & + !$acc mleVelocityMeridional, & !$acc bulkRichardsonNumberShear, & !$acc indexSurfaceVelocityZonal, & !$acc normalVelocitySurfaceLayer, & @@ -965,9 +955,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc normRelVortEdge, & !$acc surfaceVelocity, & !$acc vertVelocityTop, & - !$acc eddyVelocityX, & - !$acc eddyVelocityY, & - !$acc eddyVelocityZ, & + !$acc eddyVelocityX, & + !$acc eddyVelocityY, & + !$acc eddyVelocityZ, & !$acc vertNonLocalFlux, & !$acc vertViscTopOfEdge, & !$acc gradSSHMeridional, & @@ -987,7 +977,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc BruntVaisalaFreqTop, & !$acc vertAleTransportTop, & !$acc GMBolusVelocityZonal, & - !$acc mleVelocityZonal, & + !$acc mleVelocityZonal, & !$acc bulkRichardsonNumber, & !$acc gmStreamFuncTopOfEdge, & !$acc indexSSHGradientZonal, & @@ -1007,7 +997,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc vertTransportVelocityTop, & !$acc bulkRichardsonNumberBuoy, & !$acc tracersSurfaceLayerValue, & - !$acc mleVelocityMeridional, & + !$acc mleVelocityMeridional, & !$acc GMBolusVelocityMeridional, & !$acc bulkRichardsonNumberShear, & !$acc indexSurfaceVelocityZonal, & @@ -1171,9 +1161,9 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ normRelVortEdge, & surfaceVelocity, & vertVelocityTop, & - eddyVelocityX, & - eddyVelocityY, & - eddyVelocityZ, & + eddyVelocityX, & + eddyVelocityY, & + eddyVelocityZ, & vertNonLocalFlux, & vertViscTopOfEdge, & gradSSHMeridional, & @@ -1193,7 +1183,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ BruntVaisalaFreqTop, & vertAleTransportTop, & GMBolusVelocityZonal, & - mleVelocityZonal, & + mleVelocityZonal, & bulkRichardsonNumber, & gmStreamFuncTopOfEdge, & indexSSHGradientZonal, & @@ -1213,7 +1203,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ vertTransportVelocityTop, & bulkRichardsonNumberBuoy, & tracersSurfaceLayerValue, & - mleVelocityMeridional, & + mleVelocityMeridional, & GMBolusVelocityMeridional, & bulkRichardsonNumberShear, & indexSurfaceVelocityZonal, & From 7c984dd9d9fac55bfe073fe4c624357649d029bc Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 11 Aug 2022 12:37:06 -0500 Subject: [PATCH 14/28] fixes compile time errors --- components/mpas-ocean/src/Registry.xml | 4 ++-- .../src/analysis_members/mpas_ocn_moc_streamfunction.F | 5 +++-- components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F | 4 ++-- components/mpas-ocean/src/shared/mpas_ocn_gm.F | 5 +++-- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index d3143ba56fa7..f76150f18040 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2111,8 +2111,8 @@ - - + + diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F index 5a5f3f538329..c101e43f61fa 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F @@ -327,7 +327,7 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ real (kind=RKIND), dimension(:), pointer :: latCell real (kind=RKIND), dimension(:), pointer :: areaCell, binBoundaryMocStreamfunction real (kind=RKIND), dimension(:,:), pointer :: mocStreamvalLatAndDepth, mocStreamValLatAndDepthTotal - real (kind=RKIND), dimension(:,:), pointer :: mocStreamvalLatAndDepthGM + real (kind=RKIND), dimension(:,:), pointer :: mocStreamvalLatAndDepthGM, mocStreamvalLatAndDepthMLE real (kind=RKIND), dimension(:,:), pointer :: sumVertBinVelocity !!!! TRANSECT VARIABLES !!!! @@ -352,7 +352,8 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ integer :: currentRegion, i real (kind=RKIND), dimension(:,:,:), pointer :: mocStreamValLatAndDepthRegionLocal, & mocStreamvalLatAndDepthRegion, mocStreamValLatAndDepthRegionTotal, & - sumVertBinVelocityRegion, mocStreamvalLatAndDepthRegionGM + sumVertBinVelocityRegion, mocStreamvalLatAndDepthRegionGM, & + mocStreamvalLatAndDepthRegionMLE integer, dimension(:, :), pointer :: regionCellMasks, regionVertexMasks, regionsInGroup character (len=STRKIND), dimension(:), pointer :: regionNames, regionGroupNames integer, dimension(:), pointer :: nRegionsInGroup diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index b65d272d6cf7..44c38330b4ed 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -1676,7 +1676,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & #ifdef MPAS_OPENACC !$acc parallel loop gang vector & !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & - !$acc areaCell, layerThicknessEdgeFlux, normalMLEBolusVelocity, & + !$acc areaCell, layerThicknessEdgeFlux, normalMLEVelocity, & !$acc vertMLEBolusVelocityTop, invAreaCell, minLevelCell) & !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & !$acc dvEdge_temp, k) @@ -1698,7 +1698,7 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & ! Compute vertical velocity from the horizontal GM Bolus velocity div_huMLEBolus(k) = div_huMLEBolus(k) & - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & - * normalMLEBolusVelocity(k, iEdge) * invAreaCell1 + * normalMLEVelocity(k, iEdge) * invAreaCell1 end do end do ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index 6e968c908fda..9a6dd0e7d2c6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -144,7 +144,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND) :: c_mode1, resolution, dc_Lr ! Dimensions integer :: nCells, nEdges - integer, pointer :: nVertLevels + integer, pointer :: nVertLevels, indMLDedge integer, dimension(:), pointer :: nCellsArray, nEdgesArray type(mpas_pool_type), pointer :: tracersPool @@ -350,7 +350,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & slopeTaperDown = 1.0_RKIND + slopeTaperFactor*(slopeTaperDown - 1.0_RKIND) sfcTaper = min(RediKappaSfcTaper(k, cell1), RediKappaSfcTaper(k, cell2)) - if(k < indMLDedge(iEdge)) then + indMLDedge = max(indMLD(cell1),indMLD(cell2)) + if(k < indMLDedge) then sfcTaper = 0.0_RKIND else sfcTaper = 1.0_RKIND From 92334c161348e29c42276b01a57767bf7d341d6d Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 11 Aug 2022 12:46:04 -0500 Subject: [PATCH 15/28] Minor Formatting Fixes --- components/mpas-ocean/src/Registry.xml | 2 +- .../src/analysis_members/Registry_moc_streamfunction.xml | 2 +- .../src/analysis_members/mpas_ocn_moc_streamfunction.F | 2 +- components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index f76150f18040..cd5f8fc86fd5 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2865,7 +2865,7 @@ description="count of times redi limiter is invoked on a timestep" packages="gm" /> - diff --git a/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml b/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml index a2b32ebe6ed0..5acaf30acd30 100644 --- a/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml +++ b/components/mpas-ocean/src/analysis_members/Registry_moc_streamfunction.xml @@ -67,7 +67,7 @@ - Date: Thu, 11 Aug 2022 12:53:30 -0500 Subject: [PATCH 16/28] Further formatting fixes --- components/mpas-ocean/src/Registry.xml | 20 +++++++++---------- .../Registry_moc_streamfunction.xml | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index cd5f8fc86fd5..2129efa8d04d 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2840,7 +2840,7 @@ /> - - Date: Thu, 18 Aug 2022 00:18:10 -0500 Subject: [PATCH 17/28] changes for BFB --- components/mpas-ocean/driver/ocn_comp_mct.F | 2 +- components/mpas-ocean/src/Registry.xml | 14 +- .../mpas_ocn_high_frequency_output.F | 95 ++--- .../mpas_ocn_mixed_layer_depths.F | 4 + .../src/mode_forward/mpas_ocn_forward_mode.F | 6 +- .../mpas_ocn_time_integration_rk4.F | 129 ++++--- .../mpas_ocn_time_integration_si.F | 147 +++++--- .../mpas_ocn_time_integration_split.F | 158 +++++--- .../src/shared/mpas_ocn_diagnostics.F | 340 ++++++++++++------ .../mpas_ocn_eddy_parameterization_helpers.F | 243 ++++++++++--- .../mpas-ocean/src/shared/mpas_ocn_gm.F | 56 ++- .../src/shared/mpas_ocn_submesoscale_eddies.F | 73 +++- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 26 +- 13 files changed, 895 insertions(+), 398 deletions(-) diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index 2b9eb424a7e2..18b1cbe06246 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -973,7 +973,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr) - call ocn_eddy_compute_mixed_layer_depth(statePool) + call ocn_eddy_compute_mixed_layer_depth(statePool, forcingPool) if (config_use_GM .or. config_submesoscale_enable) then call ocn_eddy_compute_buoyancy_gradient() end if diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 2129efa8d04d..f93aae4fb5da 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -443,14 +443,22 @@ /> - - + /> + + \brief Computes vertical velocity from GM bolus velocity +!> \author Luke Van Roekel +!> \date August 2022 +!> \details +!> This routine computes the diagnostic vertical velocity from +!> the normal velocity due to GM bolus velocity +! +!----------------------------------------------------------------------- + subroutine ocn_diagnostic_solve_GMvel(layerThickEdgeFlux, normalGMBolusVelocity, & + vertGMBolusVelocityTop)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickEdgeFlux, & + normalGMBolusVelocity + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(out) :: & + vertGMBolusVelocityTop + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: nCells, nEdges + integer :: iCell, iVertex, iEdge + integer :: i, j, k + integer :: edgeSignOnCell_temp, eoe + real(kind=RKIND) :: invAreaCell1 + real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp + + real (kind=RKIND), dimension(:), allocatable:: div_huGMBolus + +#ifdef MPAS_OPENACC + !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus, div_huMLEBolus) +#endif + + ! Need owned cells for relativeVorticityCell + nCells = nCellsOwned + + if (.not. config_use_gm) return + + allocate(div_huGMBolus(nVertLevels)) +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & + !$acc areaCell, layerThicknessEdgeFlux, normalGMBolusVelocity, & + !$acc vertGMBolusVelocityTop, invAreaCell, minLevelCell) & + !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & + !$acc dvEdge_temp, k) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, iEdge, invAreaCell1, edgeSignOnCell_temp, dcEdge_temp, dvEdge_temp, & + !$omp k, div_huGMBolus) +#endif + do iCell = 1, nCells + div_huGMBolus(:) = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + edgeSignOnCell_temp = edgeSignOnCell(i, iCell) + dcEdge_temp = dcEdge(iEdge) + dvEdge_temp = dvEdge(iEdge) + do k = minLevelCell(iCell), maxLevelCell(iCell) + ! Compute vertical velocity from the horizontal GM Bolus velocity + div_huGMBolus(k) = div_huGMBolus(k) & + - layerThickEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & + * normalGMBolusVelocity(k, iEdge) * invAreaCell1 + end do + end do + ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. + vertGMBolusVelocityTop(1:minLevelCell(iCell)-1, iCell) = 0.0_RKIND + vertGMBolusVelocityTop(maxLevelCell(iCell)+1, iCell) = 0.0_RKIND + do k = maxLevelCell(iCell), 1, -1 + vertGMBolusVelocityTop(k,iCell) = vertGMBolusVelocityTop(k+1,iCell) - div_huGMBolus(k) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + deallocate(div_huGMBolus) + + end subroutine ocn_diagnostic_solve_GMvel!}}} + +!*********************************************************************** +! +! routine ocn_diagnostic_solve_MLEvel +! +!> \brief Computes vertical velocity from MLE bolus velocity +!> \author Luke Van Roekel +!> \date August 2022 +!> \details +!> This routine computes the diagnostic vertical velocity from +!> the normal velocity due to MLE (submesoscale eddy) bolus velocity +! +!----------------------------------------------------------------------- + subroutine ocn_diagnostic_solve_MLEvel(layerThickEdgeFlux, normalMLEVelocity, & + vertMLEBolusVelocityTop)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickEdgeFlux, & + normalMLEVelocity + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(out) :: & + vertMLEBolusVelocityTop + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: nCells, nEdges + integer :: iCell, iVertex, iEdge + integer :: i, j, k + integer :: edgeSignOnCell_temp, eoe + real(kind=RKIND) :: invAreaCell1 + real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp + + real (kind=RKIND), dimension(:), allocatable:: div_huMLEBolus + +#ifdef MPAS_OPENACC + !$acc declare device_resident(div_huMLEBolus) +#endif + + nCells = nCellsOwned + if (.not. config_submesoscale_enable) return + + allocate(div_huMLEBolus(nVertLevels)) +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & + !$acc areaCell, layerThicknessEdgeFlux, normalMLEVelocity, & + !$acc vertMLEBolusVelocityTop, invAreaCell, minLevelCell) & + !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & + !$acc dvEdge_temp, k) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, iEdge, invAreaCell1, edgeSignOnCell_temp, dcEdge_temp, dvEdge_temp, & + !$omp k, div_huMLEBolus) +#endif + do iCell = 1, nCells + div_huMLEBolus(:) = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + edgeSignOnCell_temp = edgeSignOnCell(i, iCell) + dcEdge_temp = dcEdge(iEdge) + dvEdge_temp = dvEdge(iEdge) + do k = minLevelCell(iCell), maxLevelCell(iCell) + ! Compute vertical velocity from the horizontal GM Bolus velocity + div_huMLEBolus(k) = div_huMLEBolus(k) & + - layerThickEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & + * normalMLEVelocity(k, iEdge) * invAreaCell1 + end do + end do + ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. + vertMLEBolusVelocityTop(1:minLevelCell(iCell)-1, iCell) = 0.0_RKIND + vertMLEBolusVelocityTop(maxLevelCell(iCell)+1, iCell) = 0.0_RKIND + do k = maxLevelCell(iCell), 1, -1 + vertMLEBolusVelocityTop(k,iCell) = vertMLEBolusVelocityTop(k+1,iCell) - div_huMLEBolus(k) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + deallocate(div_huMLEBolus) + + end subroutine ocn_diagnostic_solve_MLEvel!}}} + !*********************************************************************** ! ! routine ocn_diagnostic_solve_vortVel @@ -1563,8 +1772,7 @@ end subroutine ocn_diagnostic_solve_surfaceLayer !}}} !----------------------------------------------------------------------- subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & layerThicknessEdgeFlux, normalVelocity, normalTransportVelocity, & - normalGMBolusVelocity, normalMLEVelocity, vertVelocityTop, vertTransportVelocityTop, & - vertGMBolusVelocityTop, vertMLEBolusVelocityTop, relativeVorticityCell, & + vertVelocityTop, vertTransportVelocityTop, relativeVorticityCell, & divergence, kineticEnergyCell, tangentialVelocity)!{{{ implicit none @@ -1583,10 +1791,6 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & normalVelocity real (kind=RKIND), dimension(:,:), intent(in) :: & normalTransportVelocity - real (kind=RKIND), dimension(:,:), intent(in) :: & - normalGMBolusVelocity - real (kind=RKIND), dimension(:,:), intent(in) :: & - normalMLEVelocity !----------------------------------------------------------------- ! @@ -1598,10 +1802,6 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & vertVelocityTop real (kind=RKIND), dimension(:,:), intent(out) :: & vertTransportVelocityTop - real (kind=RKIND), dimension(:,:), intent(out) :: & - vertGMBolusVelocityTop - real (kind=RKIND), dimension(:,:), intent(out) :: & - vertMLEBolusVelocityTop real (kind=RKIND), dimension(:,:), intent(out) :: & relativeVorticityCell real (kind=RKIND), dimension(:,:), intent(out) :: & @@ -1624,11 +1824,10 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & real(kind=RKIND) :: invAreaCell1 real(kind=RKIND) :: dcEdge_temp, dvEdge_temp, r_tmp, weightsOnEdge_temp - real (kind=RKIND), dimension(:), allocatable:: div_hu,div_huTransport,div_huGMBolus - real (kind=RKIND), dimension(:), allocatable:: div_huMLEBolus + real (kind=RKIND), dimension(:), allocatable:: div_hu,div_huTransport #ifdef MPAS_OPENACC - !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus, div_huMLEBolus) + !$acc declare device_resident(div_hu, div_huTransport) #endif ! Need owned cells for relativeVorticityCell @@ -1671,94 +1870,6 @@ subroutine ocn_diagnostic_solve_vortVel(relativeVorticity, & ! allocate(div_hu(nVertLevels),div_huTransport(nVertLevels)) - if (config_submesoscale_enable) then - allocate(div_huMLEBolus(nVertLevels)) -#ifdef MPAS_OPENACC - !$acc parallel loop gang vector & - !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & - !$acc areaCell, layerThicknessEdgeFlux, normalMLEVelocity, & - !$acc vertMLEBolusVelocityTop, invAreaCell, minLevelCell) & - !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & - !$acc dvEdge_temp, k) -#else - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(i, iEdge, invAreaCell1, edgeSignOnCell_temp, dcEdge_temp, dvEdge_temp, & - !$omp k, div_huMLEBolus) -#endif - do iCell = 1, nCells - div_huGMBolus(:) = 0.0_RKIND - invAreaCell1 = invAreaCell(iCell) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - edgeSignOnCell_temp = edgeSignOnCell(i, iCell) - dcEdge_temp = dcEdge(iEdge) - dvEdge_temp = dvEdge(iEdge) - do k = minLevelCell(iCell), maxLevelCell(iCell) - ! Compute vertical velocity from the horizontal GM Bolus velocity - div_huMLEBolus(k) = div_huMLEBolus(k) & - - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & - * normalMLEVelocity(k, iEdge) * invAreaCell1 - end do - end do - ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. - vertMLEBolusVelocityTop(1:minLevelCell(iCell)-1, iCell) = 0.0_RKIND - vertMLEBolusVelocityTop(maxLevelCell(iCell)+1, iCell) = 0.0_RKIND - do k = maxLevelCell(iCell), 1, -1 - vertMLEBolusVelocityTop(k,iCell) = vertMLEBolusVelocityTop(k+1,iCell) - div_huGMBolus(k) - end do - end do -#ifndef MPAS_OPENACC - !$omp end do - !$omp end parallel -#endif - deallocate(div_huMLEBolus) - end if !config_submesoscale_enable - - if (config_use_gm) then - allocate(div_huGMBolus(nVertLevels)) -#ifdef MPAS_OPENACC - !$acc parallel loop gang vector & - !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & - !$acc areaCell, layerThicknessEdgeFlux, normalGMBolusVelocity, & - !$acc vertGMBolusVelocityTop, invAreaCell, minLevelCell) & - !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & - !$acc dvEdge_temp, k) -#else - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(i, iEdge, invAreaCell1, edgeSignOnCell_temp, dcEdge_temp, dvEdge_temp, & - !$omp k, div_huGMBolus) -#endif - do iCell = 1, nCells - div_huGMBolus(:) = 0.0_RKIND - invAreaCell1 = invAreaCell(iCell) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - edgeSignOnCell_temp = edgeSignOnCell(i, iCell) - dcEdge_temp = dcEdge(iEdge) - dvEdge_temp = dvEdge(iEdge) - do k = minLevelCell(iCell), maxLevelCell(iCell) - ! Compute vertical velocity from the horizontal GM Bolus velocity - div_huGMBolus(k) = div_huGMBolus(k) & - - layerThicknessEdgeFlux(k, iEdge) * edgeSignOnCell_temp * dvEdge_temp & - * normalGMBolusVelocity(k, iEdge) * invAreaCell1 - end do - end do - ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above. - vertGMBolusVelocityTop(1:minLevelCell(iCell)-1, iCell) = 0.0_RKIND - vertGMBolusVelocityTop(maxLevelCell(iCell)+1, iCell) = 0.0_RKIND - do k = maxLevelCell(iCell), 1, -1 - vertGMBolusVelocityTop(k,iCell) = vertGMBolusVelocityTop(k+1,iCell) - div_huGMBolus(k) - end do - end do -#ifndef MPAS_OPENACC - !$omp end do - !$omp end parallel -#endif - deallocate(div_huGMBolus) - end if !config_use_gm - #ifdef MPAS_OPENACC !$acc parallel loop gang vector & !$acc present(divergence, kineticEnergyCell, & @@ -3194,7 +3305,8 @@ subroutine ocn_reconstruct_eddy_vectors(meshPool) !{{{ transportVelocityMeridional & ) - call mpas_reconstruct(meshPool, normalGMBolusVelocity, & + if(config_use_gm) then + call mpas_reconstruct(meshPool, normalGMBolusVelocity, & eddyVelocityX, & eddyVelocityY, & eddyVelocityZ, & @@ -3202,7 +3314,17 @@ subroutine ocn_reconstruct_eddy_vectors(meshPool) !{{{ GMBolusVelocityMeridional & ) - call mpas_reconstruct(meshPool, normalMLEvelocity, & + call mpas_reconstruct(meshPool, gmStreamFuncTopOfEdge, & + GMStreamFuncX, & + GMStreamFuncY, & + GMStreamFuncZ, & + GMStreamFuncZonal, & + GMStreamFuncMeridional & + ) + end if + + if (config_submesoscale_enable) then + call mpas_reconstruct(meshPool, normalMLEvelocity, & eddyVelocityX, & eddyVelocityY, & eddyVelocityZ, & @@ -3210,13 +3332,7 @@ subroutine ocn_reconstruct_eddy_vectors(meshPool) !{{{ mleVelocityMeridional & ) - call mpas_reconstruct(meshPool, gmStreamFuncTopOfEdge, & - GMStreamFuncX, & - GMStreamFuncY, & - GMStreamFuncZ, & - GMStreamFuncZonal, & - GMStreamFuncMeridional & - ) + end if call mpas_timer_stop('reconstruct eddy vecs') diff --git a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F index d8edd70eaa5e..fb55d6b5df24 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F @@ -166,8 +166,8 @@ subroutine ocn_eddy_compute_buoyancy_gradient() do iEdge = 1, nEdges if (maxLevelEdgeTop(iEdge) .GE. minLevelEdgeBot(iEdge)) then do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1 - gradBuoyEddy(k,iEdge) = gravity*(gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & - * gradZMidTopOfEdge(k,iEdge))/rho_sw + gradBuoyEddy(k,iEdge) = (gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & + * gradZMidTopOfEdge(k,iEdge)) end do end if end do @@ -193,79 +193,218 @@ end subroutine ocn_eddy_compute_buoyancy_gradient !> The calculation is moved out of analysis members for this reason !*********************************************************************** - subroutine ocn_eddy_compute_mixed_layer_depth(statePool) + subroutine ocn_eddy_compute_mixed_layer_depth(statePool, forcingPool) - type (mpas_pool_type), pointer, intent(in) :: statePool + type (mpas_pool_type), pointer, intent(in) :: statePool, forcingPool integer :: iEdge, nEdges, refIndex, cell1, cell2, nCells, k, iCell real(kind=RKIND) :: dDenThres, den_ref_lev real(kind=RKIND),dimension(:), allocatable :: depth integer, dimension(:), pointer :: landIceMask real (kind=RKIND), dimension(:,:), pointer :: layerThickness + real (kind=RKIND), dimension(:), pointer :: landIceDraft, landIcePressure logical :: found_den_mld - real (kind=RKIND) :: dV, dVp1, refDepth, coeffs(2), mldTemp + real (kind=RKIND) :: dV, dVp1, refDepth, coeffs(2), mldTemp, dz + real (kind=RKIND), dimension(:,:), allocatable :: pressureAdjustedForLandIce - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) nCells = nCellsAll - allocate(depth(nVertLevels)) - refDepth = config_mld_reference_depth + if ( config_eddyMLD_use_old ) then + call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) + call mpas_pool_get_array(forcingPool, 'landIceDraft', landIceDraft) + call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask) + + allocate(pressureAdjustedForLandIce(nVertLevels, size(pressure,2))) + + pressureAdjustedForLandIce(:,:) = pressure(:,:) + !If landice cavity remove land ice pressure to search for ML depth + if ( associated(landIcePressure) ) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1,nCells + do k = 1,maxLevelCell(iCell) + pressureAdjustedForLandIce(k,iCell) = pressureAdjustedForLandIce(k,iCell) & + - landIcePressure(iCell) + end do + end do + !$omp end do + !$omp end parallel + end if - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(depth, refIndex, coeffs, found_den_mld, k, den_ref_lev, dVp1, dV, mldTemp) - do iCell = 1,nCells + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(refIndex, found_den_mld, k, localvals, den_ref_lev, dVp1, dV, mldTemp) + do iCell = 1,nCells - found_den_mld = .false. + !Initialize RefIndex for cases of very shallow columns + refIndex = maxLevelCell(iCell) + found_den_mld = .false. - depth(1) = 0.5_RKIND*layerThickness(1,iCell) - do k=2,maxLevelCell(iCell) - depth(k) = depth(k-1) + 0.5_RKIND*(layerThickness(k-1,iCell) + & - layerThickness(k ,iCell)) - end do + do k=1, maxLevelCell(iCell)-1 + if(pressureAdjustedForLandIce(k+1,iCell) > config_eddyMLD_reference_pressure) then - do k=1, maxLevelCell(iCell)-1 - if(depth(k+1) >= refDepth) then - coeffs(2) = (potentialDensity(k+1,iCell) - potentialDensity(k,iCell)) & - / (depth(k+1) - depth(k)) - coeffs(1) = potentialDensity(k,iCell) - coeffs(2)*depth(k) + coeffs(2) = (potentialDensity(k+1,iCell) - potentialDensity(k,iCell)) / & + (pressureAdjustedForLandIce(k+1,iCell) - pressureAdjustedForLandIce(k,iCell)) + coeffs(1) = potentialDensity(k,iCell) - coeffs(2)*pressureAdjustedForLandIce(k,iCell) - den_ref_lev = coeffs(2)*refDepth + coeffs(1) + den_ref_lev = coeffs(2)*config_eddyMLD_reference_pressure + coeffs(1) - refIndex=k - exit - end if - end do + refIndex=k + exit + end if + end do + + do k=refIndex,maxLevelCell(iCell)-1 + + if(.not. found_den_mld .and. abs(potentialDensity(k+1,iCell) - den_ref_lev) .ge. & + config_eddyMLD_dens_threshold) then + dVp1 = abs(potentialDensity(k+1,iCell) - den_ref_lev) + dV = abs(potentialDensity(k ,iCell) - den_ref_lev) + + coeffs(2) = (zMid(k+1,iCell) - zMid(k,iCell)) / (dVp1 - dV) + coeffs(1) = zMid(k,iCell) - coeffs(2)*dV - do k=refIndex,maxLevelCell(iCell)-1 + mldTemp = coeffs(2)*config_eddyMLD_dens_threshold + coeffs(1) - if(.not. found_den_mld .and. (potentialDensity(k+1,iCell) - den_ref_lev) .ge. & - config_mixedLayerDepths_crit_dens_threshold) then - dVp1 = (potentialDensity(k+1,iCell) - den_ref_lev) - dV = (potentialDensity(k ,iCell) - den_ref_lev) + mldTemp=max(mldTemp,zMid(k+1,iCell)) !make sure MLD isn't deeper than zMid(k+1) + dThreshMLD(iCell)=abs(min(mldTemp,zMid(k,iCell))) !MLD should be deeper than zMid(k) + indMLD(iCell) = k+1 + found_den_mld = .true. + exit + end if + end do - coeffs(2) = (depth(k+1) - depth(k)) / (dVp1 - dV) - coeffs(1) = depth(k) - coeffs(2)*dV + ! if the mixed layer depth is not found, it is set to the depth of the bottom most level - mldTemp = coeffs(2)*config_mixedLayerDepths_crit_dens_threshold + coeffs(1) - mldTemp=min(mldTemp,depth(k+1)) !make sure MLD isn't deeper than zMid(k+1) - dThreshMLD(iCell)=abs(max(mldTemp,depth(k))) !MLD should be deeper than zMid(k) - indMLD(iCell) = k+1 - found_den_mld = .true. - exit + if(.not. found_den_mld) then + dThreshMLD(iCell) = abs(zMid(maxLevelCell(iCell),iCell)) + indMLD(iCell) = maxLevelCell(iCell) end if - end do + end do !iCell + !$omp end do + !$omp end parallel + + if (associated(landIceMask) ) then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + dThreshMLD(iCell) = dThreshMLD(iCell) - abs(landIceDraft(iCell))*landIceMask(iCell) + end do + !$omp end do + !$omp end parallel + end if - ! if the mixed layer depth is not found, it is set to the depth of the bottom most level - if(.not. found_den_mld) then - dThreshMLD(iCell) = depth(maxLevelCell(iCell)) - indMLD(iCell) = maxLevelCell(iCell) - end if - end do !iCell - !$omp end do - !$omp end parallel + deallocate(pressureAdjustedForLandIce) + + else ! fixed mixed layer depth code + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) + allocate(depth(nVertLevels)) + refDepth = config_eddyMLD_reference_depth + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(depth, refIndex, coeffs, found_den_mld, k, den_ref_lev, dVp1, dV, mldTemp) + do iCell = 1,nCells + + found_den_mld = .false. + + depth(1) = 0.5_RKIND*layerThickness(1,iCell) + do k=2,maxLevelCell(iCell) + depth(k) = depth(k-1) + 0.5_RKIND*(layerThickness(k-1,iCell) + & + layerThickness(k ,iCell)) + end do + + do k=1, maxLevelCell(iCell)-1 + if(depth(k+1) >= refDepth) then + coeffs(2) = (potentialDensity(k+1,iCell) - potentialDensity(k,iCell)) & + / (depth(k+1) - depth(k)) + coeffs(1) = potentialDensity(k,iCell) - coeffs(2)*depth(k) + + den_ref_lev = coeffs(2)*refDepth + coeffs(1) + + refIndex=k + exit + end if + end do + + do k=refIndex,maxLevelCell(iCell)-1 + + if(.not. found_den_mld .and. (potentialDensity(k+1,iCell) - den_ref_lev) .ge. & + config_eddyMLD_dens_threshold) then + dVp1 = (potentialDensity(k+1,iCell) - den_ref_lev) + dV = (potentialDensity(k ,iCell) - den_ref_lev) + + coeffs(2) = (depth(k+1) - depth(k)) / (dVp1 - dV) + coeffs(1) = depth(k) - coeffs(2)*dV + + mldTemp = coeffs(2)*config_eddyMLD_dens_threshold + coeffs(1) + mldTemp=min(mldTemp,depth(k+1)) !make sure MLD isn't deeper than zMid(k+1) + dThreshMLD(iCell)=abs(max(mldTemp,depth(k))) !MLD should be deeper than zMid(k) + indMLD(iCell) = k+1 + found_den_mld = .true. + exit + end if + end do + + ! if the mixed layer depth is not found, it is set to the depth of the bottom most level + if(.not. found_den_mld) then + dThreshMLD(iCell) = depth(maxLevelCell(iCell)) + indMLD(iCell) = maxLevelCell(iCell) + end if + end do !iCell + !$omp end do + !$omp end parallel + + deallocate(depth) + + end if !end mld calculation choice + +end subroutine ocn_eddy_compute_mixed_layer_depth + +!*********************************************************************** +! +! routine interp_mld +! +!> \brief Interpolates between model layers +!> \author Luke Van Roekel +!> \date September 2015 +!> \details +!> This routine conducts computations to compute various field values +!> between model levels (in pressure or depth) or could interpolate +!> between temperature/salinity/density values. Interpolations are +!> of the form +!> y = coeffs(1)*x^3 + coeffs(2)*x^2 + coeffs(3)*x + coeffs(4) +! +!----------------------------------------------------------------------- + + subroutine interp_mld(y0,y1,x0,x1,xT,interp_f,yT,xm1,ym1) + + integer,intent(in) :: interp_f ! linear, quadratic, or spline + real(kind=RKIND),intent(in) :: y0,y1,x0,x1,xT + real(kind=RKIND),intent(inout) :: yT + real(kind=RKIND),optional,intent(in) :: xm1,ym1 + ! these values are to match the slope at a given point + +!------------------------------------------------------------------------ +! +! Local variables for the interpolations +! +!------------------------------------------------------------------------ + + real(kind=RKIND) :: coeffs(4) ! stores the coefficients for the interp + real(kind=RKIND) :: Minv(4,4) ! holds values for computing quad and spline + real(kind=RKIND) :: det + real(kind=RKIND) :: rhs(4) + integer :: k,k2 + + coeffs(:) = 0.0_RKIND + Minv(:,:) = 0.0_RKIND + rhs(:) = 0.0_RKIND + + coeffs(2) = (y1-y0)/(x1-x0) + coeffs(1) = y0 - coeffs(2)*x0 - deallocate(depth) - end subroutine ocn_eddy_compute_mixed_layer_depth + yT = coeffs(4)*xT**3 + coeffs(3)*xT**2 + coeffs(2)*xT + coeffs(1) + end subroutine interp_mld end module ocn_eddy_parameterization_helpers diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index 9a6dd0e7d2c6..905298d5a092 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -35,7 +35,8 @@ module ocn_gm !-------------------------------------------------------------------- public :: ocn_GM_compute_Bolus_velocity, & - ocn_GM_init + ocn_GM_init, & + ocn_GM_add_to_transport_vel !-------------------------------------------------------------------- ! @@ -144,7 +145,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND) :: c_mode1, resolution, dc_Lr ! Dimensions integer :: nCells, nEdges - integer, pointer :: nVertLevels, indMLDedge + integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray, nEdgesArray type(mpas_pool_type), pointer :: tracersPool @@ -350,12 +351,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & slopeTaperDown = 1.0_RKIND + slopeTaperFactor*(slopeTaperDown - 1.0_RKIND) sfcTaper = min(RediKappaSfcTaper(k, cell1), RediKappaSfcTaper(k, cell2)) - indMLDedge = max(indMLD(cell1),indMLD(cell2)) - if(k < indMLDedge) then - sfcTaper = 0.0_RKIND - else - sfcTaper = 1.0_RKIND - end if sfcTaperUp = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND) sfcTaperDown = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND) @@ -715,7 +710,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* & + gradBuoyEddy(k,iEdge) ! Second to next to the last rows do k = minLevelEdgeBot(iEdge) + 2, maxLevelEdgeTop(iEdge) - 1 @@ -727,7 +723,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* & + gradBuoyEddy(k,iEdge) end do ! Last row @@ -738,7 +735,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge)) tridiagB(k - 1) = -2.0_RKIND*cGMphaseSpeed(iEdge)**2/(layerThickEdgeMean(k - 1, iEdge) & *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge - rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gradBuoyEddy(k,iEdge) + rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* & + gradBuoyEddy(k,iEdge) ! Total number of rows N = maxLevelEdgeTop(iEdge) - minLevelEdgeBot(iEdge) @@ -838,6 +836,42 @@ subroutine tridiagonal_solve(a, b, c, r, x, n) !{{{ end subroutine tridiagonal_solve !}}} +!*********************************************************************** +! +! routine ocn_GM_add_to_transport_velocity +! +!> \brief Adds GM to transportVel +!> \details +!> Adds the GM bolus velocity to the normal transport velocity if enabled +! +!----------------------------------------------------------------------- + + subroutine ocn_GM_add_to_transport_vel(normalTransportVelocity, & + nEdges, nVertLevels)!{{{ + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + normalTransportVelocity + + integer, intent(in) :: nEdges, nVertLevels + + integer :: iEdge, k + + if (.not. config_use_gm) return + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iEdge = 1, nEdges + do k = 1, nVertLevels + normalTransportVelocity(k,iEdge) = normalTransportVelocity(k,iEdge) + & + normalGMBolusVelocity(k,iEdge) + end do + end do + !$omp end do + !$omp end parallel + + end subroutine ocn_GM_add_to_transport_vel!}}} + !*********************************************************************** ! ! routine ocn_GM_init diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F index 1dbf51b94fd4..e2e5737127bb 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F @@ -36,10 +36,12 @@ module ocn_submesoscale_eddies !-------------------------------------------------------------------- public :: ocn_submesoscale_compute_velocity, & - ocn_submesoscale_init + ocn_submesoscale_init, & + ocn_submesoscale_finalize, & + ocn_MLE_add_to_transport_vel real(kind=RKIND) :: tau, LfMin, Ce, dsMax - + real(kind=RKIND),dimension(:),allocatable :: time_scale !*********************************************************************** @@ -90,6 +92,10 @@ subroutine ocn_submesoscale_compute_velocity() allocate(streamFunction(nVertLevels+1)) allocate(zEdge(nVertLevels+1)) + !compute the mixed layer average buoyancy gradient and bvf**2 + !This is thickness weighted to allow for non uniform grid spacing + !NOTE, it is assumed that the value of bvf**2 at the bottom of cell 1 + ! applies to the entire first cell !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1,cell2,bvfML,zEdge,streamFunction,gradBuoyML,hML,bvfAv, & @@ -102,6 +108,7 @@ subroutine ocn_submesoscale_compute_velocity() streamFunction(:) = 0.0_RKIND bvfML = 0.0_RKIND gradBuoyML = 0.0_RKIND + !the BFV is defined at cell interfaces in the vertical so there is an offset in index bvfAv = sqrt(0.5_RKIND*(max(BruntVaisalaFreqTop(minLevelCell(cell1)+1,cell1),1.0E-20_RKIND) + & max(BruntVaisalaFreqTop(minLevelCell(cell2)+1,cell2),1.0E-20_RKIND))) hML = 0.5_RKIND*layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge) @@ -114,10 +121,11 @@ subroutine ocn_submesoscale_compute_velocity() end if indMLDedge = min(indMLD(cell1),indMLD(cell2)) - do k = minLevelEdgeBot(iEdge)+2,indMLDedge + do k = minLevelEdgeBot(iEdge)+1,indMLDedge hAv = 0.5_RKIND*(layerThickEdgeMean(k,iEdge) + layerThickEdgeMean(k-1,iEdge)) bvfML = bvfML + hAv*bvfAv hML = hML + hAv + !the bvfAV is updated after bvfML due to the offset in index from cell center versus cell top bvfAv = sqrt(0.5_RKIND*(max(BruntVaisalaFreqTop(k,cell1),1.0E-20_RKIND) + & max(BruntVaisalaFreqTop(k,cell2),1.0E-20_RKIND))) gradBuoyML = gradBuoyML + hAv*gradBuoyEddy(k,iEdge) @@ -132,15 +140,15 @@ subroutine ocn_submesoscale_compute_velocity() end do zMLD = 0.5_RKIND*(dThreshMLD(cell1)+dThreshMLD(cell2)) - Lf = max(Lfmin, abs(gradBuoyML)*zMLD / (1.0E-20_RKIND + fEdge(iEdge)**2), & - bvfML*zMLD / abs(fEdge(iEdge))) + Lf = max(Lfmin, abs(gradBuoyML)*zMLD / time_scale(iEdge)**2.0, & + bvfML*zMLD / time_scale(iEdge)) ds = min(dcEdge(iEdge),dsMax) do k = minLevelEdgeTop(iEdge)+1,maxLevelEdgeTop(iEdge) mu = max(0.0_RKIND,(1.0_RKIND - (2.0_RKIND*zEdge(k) / zMLD + 1.0_RKIND)**2.0)* & (1.0_RKIND + 5.0_RKIND/21.0_RKIND*(2.0_RKIND*zEdge(k)/zMLD + 1.0_RKIND)**2.0)) - streamFunction(k) = Ce*ds/Lf*zMLD**2.0*gradBuoyML/sqrt(fEdge(iEdge)**2.0 + tau**(-2.0))*mu + streamFunction(k) = Ce*ds/Lf*zMLD**2.0*gradBuoyML/time_scale(iEdge)*mu end do ! integrate in vertical to get the velocity @@ -157,6 +165,41 @@ subroutine ocn_submesoscale_compute_velocity() end subroutine ocn_submesoscale_compute_velocity + !*********************************************************************** + ! + ! routine ocn_add_MLE_to_transport_vel + ! + !> \brief Submesoscale parameterization add to transport Vel + !> \details + !> adds the MLE induced velocity to the transport velocity + !> + !*********************************************************************** + + subroutine ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) + + real(kind=RKIND), dimension(:,:), intent(inout) :: & + normalTransportVelocity + + integer, intent(in) :: nEdges + + integer :: iEdge, k + + if (.not. config_submesoscale_enable) return + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iEdge = 1, nEdges + do k = 1, nVertLevels + normalTransportVelocity(k,iEdge) = normalTransportVelocity(k,iEdge) + & + normalMLEvelocity(k,iEdge) + end do + end do + !$omp end do + !$omp end parallel + + end subroutine ocn_MLE_add_to_transport_vel + !*********************************************************************** ! ! routine ocn_submesoscale_init @@ -171,12 +214,30 @@ subroutine ocn_submesoscale_init(err) integer, intent(out) :: err !< Output: error flag + integer :: iEdge + err = 0 Ce = config_submesoscale_ce Lfmin = config_submesoscale_lfmin dsmax = config_submesoscale_ds_max tau = config_submesoscale_tau + allocate(time_scale(nEdgesAll)) + + !$omp parallel + !$omp do schedule(runtime) + do iEdge=1,nEdgesAll + time_scale(iEdge) = sqrt(fEdge(iEdge)**2 + tau**(-2.0)) + end do + !$omp end do + !$omp end parallel + end subroutine ocn_submesoscale_init + subroutine ocn_submesoscale_finalize() + + deallocate(time_scale) + + end subroutine ocn_submesoscale_finalize + end module ocn_submesoscale_eddies diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F index 8c70d571613d..47d17483e276 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -238,10 +238,14 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea k = minLevelEdgeBot(iEdge) kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & - RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, 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)) + 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)) do iTr = 1, nTracers tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) @@ -276,7 +280,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea 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)) + 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)) do iTr = 1, nTracers ! \kappa_2 \nabla \phi on edge @@ -313,7 +319,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea k = maxLevelEdgeTop(iEdge) kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & - RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, 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)) do iTr = 1, nTracers tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) @@ -368,8 +376,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea ! 2.0 in next line is because a dot product on a C-grid ! requires a factor of 1/2 to average to the cell center. flux_term3 = rediKappaCell(iCell)*2.0_RKIND* & - (RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) & - * invAreaCell - RediKappaScaling(k + 1, iCell) & + (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) & + * invAreaCell - RediKappaSfcTaper(k + 1, iCell)*RediKappaScaling(k + 1, iCell) & *fluxRediZTop(iTr, k + 1) * invAreaCell) redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3 @@ -438,9 +446,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* & - (RediKappaScaling(k, iCell)* & - redi_term3_topOfCell(iTr, k, iCell) - & - RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell)) + (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)* & + redi_term3_topOfCell(iTr, k, iCell) - RediKappaSfcTaper(k + 1, iCell) & + *RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell)) end do end do From 4130dd9ff0d5517e55f72d3ed50a3410febf5c94 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Thu, 18 Aug 2022 07:22:20 -0600 Subject: [PATCH 18/28] small corrections --- .../src/mode_forward/mpas_ocn_time_integration_rk4.F | 8 ++++---- .../src/mode_forward/mpas_ocn_time_integration_si.F | 4 ++-- .../src/shared/mpas_ocn_eddy_parameterization_helpers.F | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 408cc6ed5098..377a1ce338ef 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -794,7 +794,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ ! add submesoscale aend gm components if requested call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(normalTransportVelocity) + call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! ------------------------------------------------------------------ ! End: Accumulating various parameterizations of the transport velocity @@ -964,7 +964,7 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh ! add submesoscale and GM contributions if requested call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(normalTransportVelocity) + call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. if (associated(highFreqThicknessProvis)) then @@ -1036,7 +1036,7 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig !$omp end parallel call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(normalTransportVelocity) + call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. if (associated(highFreqThicknessProvis)) then @@ -1387,7 +1387,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, ! add submesoscale and GM contributions if requested call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(normalTransportVelocity) + call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! ------------------------------------------------------------------ ! End: Accumulating various parametrizations of the transport velocity diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 29784e76851f..d64287385d2b 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -1684,7 +1684,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !add GM and submesoscale contributions if requested call ocn_GM_add_to_transport_vel(uTemp, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(uTemp) + call ocn_MLE_add_to_transport_vel(uTemp, nEdges) !$omp parallel !$omp do schedule(runtime) & @@ -1701,7 +1701,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! add GM and submesoscale contributions if requested call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) - call ocn_MLE_add_to_transport_vel(normalTransportVelocity) + call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) !$omp parallel !$omp do schedule(runtime) & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F index fb55d6b5df24..236da8d38849 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F @@ -233,7 +233,7 @@ subroutine ocn_eddy_compute_mixed_layer_depth(statePool, forcingPool) !$omp parallel !$omp do schedule(runtime) & - !$omp private(refIndex, found_den_mld, k, localvals, den_ref_lev, dVp1, dV, mldTemp) + !$omp private(refIndex, found_den_mld, k, coeffs, den_ref_lev, dVp1, dV, mldTemp) do iCell = 1,nCells !Initialize RefIndex for cases of very shallow columns From d62adba73f4566816db89fb903ad204124af32dd Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 18 Aug 2022 13:29:18 -0500 Subject: [PATCH 19/28] Formatting Fixes --- components/mpas-ocean/src/Registry.xml | 2 +- components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F | 2 +- .../mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F | 9 ++++----- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index f93aae4fb5da..0d63101a29a3 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -358,7 +358,7 @@ /> Date: Sat, 20 Aug 2022 23:08:54 -0500 Subject: [PATCH 20/28] updates comments in submesoscale branch --- .../src/shared/mpas_ocn_submesoscale_eddies.F | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F index 27a1337ffafa..99bab62abb5e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F @@ -92,10 +92,13 @@ subroutine ocn_submesoscale_compute_velocity() allocate(streamFunction(nVertLevels+1)) allocate(zEdge(nVertLevels+1)) - !compute the mixed layer average buoyancy gradient and bvf**2 - !This is thickness weighted to allow for non uniform grid spacing - !NOTE, it is assumed that the value of bvf**2 at the bottom of cell 1 - ! applies to the entire first cell + !compute the mixed layer average buoyancy gradient and bvf + !This is thickness weighted to allow for non uniform grid spacing, + !and centered on layer interfaces. The mixed layer goes from the + !ocean top layer to mid-depth of layer(indMLDedge) + !NOTE, gradBuoyEddy and BFV are defined at the top cell interface + !and not valid for minLevelCell, so properties for the top layer + !rely on the values at (minLevelCell+1). !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1,cell2,bvfML,zEdge,streamFunction,gradBuoyML,hML,bvfAv, & @@ -109,6 +112,7 @@ subroutine ocn_submesoscale_compute_velocity() bvfML = 0.0_RKIND gradBuoyML = 0.0_RKIND !the BFV is defined at cell interfaces in the vertical so there is an offset in index + !this is the top half layer, assuming homogeneity of top layer bvfAv = sqrt(0.5_RKIND*(max(BruntVaisalaFreqTop(minLevelCell(cell1)+1,cell1),1.0E-20_RKIND) + & max(BruntVaisalaFreqTop(minLevelCell(cell2)+1,cell2),1.0E-20_RKIND))) hML = 0.5_RKIND*layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge) From 8b1b9ce60cb2bec1aa42b9150535dd54026d6e77 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Mon, 22 Aug 2022 08:27:37 -0600 Subject: [PATCH 21/28] Update E3SM namelist scripts --- components/mpas-ocean/bld/build-namelist | 6 ++++-- .../mpas-ocean/bld/build-namelist-section | 6 ++++-- .../namelist_defaults_mpaso.xml | 6 ++++-- .../namelist_definition_mpaso.xml | 20 +++++++++++++++++-- components/mpas-ocean/src/Registry.xml | 4 ++-- 5 files changed, 32 insertions(+), 10 deletions(-) diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 6011661e6ef3..74d6932ab8e3 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -606,8 +606,10 @@ add_default($nl, 'config_GMRedi_Rossby_ramp_max'); # Namelist group: eddy_parameterization # ######################################### -add_default($nl, 'config_mixedLayerDepths_crit_dens_threshold'); -add_default($nl, 'config_mld_reference_depth'); +add_default($nl, 'config_eddyMLD_dens_threshold'); +add_default($nl, 'config_eddyMLD_reference_depth'); +add_default($nl, 'config_eddyMLD_reference_pressure'); +add_default($nl, 'config_eddyMLD_use_old'); #################################### # Namelist group: Rayleigh_damping # diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index c1c185840e24..e0f2449e52a6 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -149,8 +149,10 @@ add_default($nl, 'config_GMRedi_Rossby_ramp_max'); # Namelist group: eddy_parameterization # ######################################### -add_default($nl, 'config_mixedLayerDepths_crit_dens_threshold'); -add_default($nl, 'config_mld_reference_depth'); +add_default($nl, 'config_eddyMLD_dens_threshold'); +add_default($nl, 'config_eddyMLD_reference_depth'); +add_default($nl, 'config_eddyMLD_reference_pressure'); +add_default($nl, 'config_eddyMLD_use_old'); #################################### # Namelist group: Rayleigh_damping # diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 6df175e82fda..731e78fa3248 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -199,8 +199,10 @@ 3.0 -0.03 -10 +0.03 +10 +1.0e5 +.true. .false. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 3943d0e501ae..6d230a7889d3 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -699,7 +699,7 @@ Default: Defined in namelist_defaults.xml - potential density change relative to surface for mixed layer depth threshold method. This calculation is used for the Redi tapering, GM N2_dependent bolus kappa, and the submesoscale eddy parameterization @@ -707,7 +707,7 @@ Valid values: suggested range 0.01 less than or equal to thresh less than or equ Default: Defined in namelist_defaults.xml - reference depth for threshold computation @@ -715,6 +715,22 @@ Valid values: any positive real, near 10 Default: Defined in namelist_defaults.xml + +reference pressure for original mixed layer depth calculation + +Valid values: positive reals around 1.0e5 +Default: Defined in namelist_defaults.xml + + + +switches from old dThreshMLD calculation to new (fixed one) + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 0d63101a29a3..ae6de2b6dfce 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -356,11 +356,11 @@ description="efficiency of submesoscale eddies" possible_values="0.06 - 0.08" /> - - From 0404c7ee03b5ecb7125a70be0eedd5a18132b95d Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 22 Aug 2022 10:21:14 -0500 Subject: [PATCH 22/28] Update bld files to match Registry changes using autogenerating scripts --- .../mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 6d230a7889d3..592dc3813320 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -529,7 +529,7 @@ Default: Defined in namelist_defaults.xml category="submesoscale_eddy_parameterization" group="submesoscale_eddy_parameterization"> minimum frontal width (meters) -Valid values: between 500 and 2000m +Valid values: between 200 and 5000m Default: Defined in namelist_defaults.xml From 37b6b94465bba2afbf8821bed2344f23cbe90ebb Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Wed, 24 Aug 2022 10:46:48 -0600 Subject: [PATCH 23/28] update outdated comments --- .../src/mode_forward/mpas_ocn_time_integration_split.F | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 6902e8143171..7f1ac1f921e8 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -1569,7 +1569,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ allocate(uTemp(nVertLevels,nEdges)) !Compute uTemp - ! velocity for normalVelocityCorrectionection is + ! velocity for normalVelocityCorrection is ! normalBarotropicVelocity + ! normalBaroclinicVelocity + uBolus !$omp parallel @@ -1633,13 +1633,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! normalTransportVelocity = normalBarotropicVelocity ! + normalBaroclinicVelocity ! + normalGMBolusVelocity + ! + normalMLEvelocity ! + normalVelocityCorrection - ! This is u used in advective terms for layerThickness + ! This is the velocity used in advective terms for layerThickness ! and tracers in tendency calls in stage 3. - !mrp note: in QC version, there is an if - ! (config_use_GM) on adding normalGMBolusVelocity - ! I think it is not needed because - ! normalGMBolusVelocity=0 when GM not on. normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge)*( & normalTransportVelocity(k,iEdge) + & normalVelocityCorrection) From fa67aa9e2d8357252e447eb8f371f748cdf14840 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Wed, 24 Aug 2022 08:05:07 -0600 Subject: [PATCH 24/28] RK4 change for bfb match intel opt --- .../mpas_ocn_time_integration_rk4.F | 32 +++++++++++-------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 377a1ce338ef..160dfce97950 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -954,13 +954,15 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) - !$omp parallel - !$omp do schedule(runtime) - do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) - end do - !$omp end do - !$omp end parallel + if (config_use_GM .or. config_submesoscale_enable) then + !$omp parallel + !$omp do schedule(runtime) + do iEdge = 1, nEdges + normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + end do + !$omp end do + !$omp end parallel + end if ! add submesoscale and GM contributions if requested call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) @@ -1027,13 +1029,15 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) - !$omp parallel - !$omp do schedule(runtime) - do iEdge = 1, nEdges - normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) - end do - !$omp end do - !$omp end parallel + if (config_use_GM .or. config_submesoscale_enable) then + !$omp parallel + !$omp do schedule(runtime) + do iEdge = 1, nEdges + normalTransportVelocity(:, iEdge) = normalVelocityProvis(:, iEdge) + end do + !$omp end do + !$omp end parallel + end if call ocn_GM_add_to_transport_vel(normalTransportVelocity, nEdges, nVertLevels) call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) From 2c6f496a2bbc764434817609c3d8110287fb828c Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Wed, 24 Aug 2022 11:01:45 -0600 Subject: [PATCH 25/28] Split Explicit change for bfb match, intel Opt. Revert later --- components/mpas-ocean/src/Registry.xml | 1 - .../mpas_ocn_time_integration_split.F | 17 +++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index ae6de2b6dfce..4125943d6b77 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -2844,7 +2844,6 @@ /> Date: Wed, 24 Aug 2022 12:08:10 -0600 Subject: [PATCH 26/28] Registry formatting fix --- components/mpas-ocean/src/Registry.xml | 34 +++++++++++++------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 4125943d6b77..33d202b1cc2f 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -330,7 +330,7 @@ description="Redi diagonal terms (2 and 3) are turned off from layer 1 through config_Redi_min_layers_diag_terms-1, and on from config_Redi_min_layers_diag_terms to nVertLevels. The Redi diagonal terms are not guaranteed to produce bounded tracer fields, and in practice produce growing temperatures in a few columns with fewer than 5 vertical cells. Redi is meant for isopycnal mixing in the deep ocean, so not applying Redi diagonal terms in very shallow regions is an acceptable solution." possible_values="any integer between 0 (all layers on) and nVertLevels (all layers off)" /> - @@ -421,7 +421,7 @@ description="factor multiplying the Rhines length in the scheme from Eden Greatbatch (2008) Ocean Modeling -- Equation (28)" possible_values="small positive values less than equal to one" /> - @@ -433,23 +433,23 @@ description="Maximum value in grid cell size for GM $\kappa$ ramp function. Here cell size refers to dcEdge. Used when config_GM_horizontal_taper is set to ramp." possible_values="Any positive real value." /> - - - - + - + - - - @@ -3028,11 +3028,11 @@ description="Redi Kappa value. On output, it has already been multiplied by the horizontal taper array RediHorizontalTaper (as opposed to gmBolusKappa, which has not been multiplied by the horizontal taper)." packages="gm" /> - - From 9e823203868e846ec3657bf2c20368a7aebda3bf Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Fri, 26 Aug 2022 10:19:00 -0600 Subject: [PATCH 27/28] Fix leftover k+1 in MOC from bottom-up calculation --- .../analysis_members/mpas_ocn_moc_streamfunction.F | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F index f61a8ddf5c06..fadd74d8f86e 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_moc_streamfunction.F @@ -510,8 +510,8 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ end do do k = 2, nVertLevels mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & - mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & - + sumTransport(k + 1, iTransect) + mocStreamValLatAndDepthRegionLocal(1, k - 1, iTransect) & + + sumTransport(k - 1, iTransect) end do end do @@ -608,8 +608,8 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ end do do k = 2, nVertLevels mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & - mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & - + sumTransport(k + 1, iTransect) + mocStreamValLatAndDepthRegionLocal(1, k - 1, iTransect) & + + sumTransport(k - 1, iTransect) end do end do @@ -714,8 +714,8 @@ subroutine ocn_compute_moc_streamfunction(domain, timeLevel, err)!{{{ end do do k = 2, nVertLevels mocStreamValLatAndDepthRegionLocal(1, k, iTransect) = & - mocStreamValLatAndDepthRegionLocal(1, k + 1, iTransect) & - + sumTransport(k + 1, iTransect) + mocStreamValLatAndDepthRegionLocal(1, k - 1, iTransect) & + + sumTransport(k - 1, iTransect) end do end do From f1b83292a534862f8653a86524071e0207751b74 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Sun, 28 Aug 2022 16:09:55 -0500 Subject: [PATCH 28/28] Updates from Az to get this PR working on OpenACC --- .../src/mode_forward/mpas_ocn_time_integration_rk4.F | 1 + components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 160dfce97950..08496ef86225 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -1320,6 +1320,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight, end if if (config_submesoscale_enable) then !$acc update device (normalMLEvelocity) + end if !$acc enter data copyin(atmosphericPressure, seaIcePressure) !$acc enter data copyin(sshCur) !$acc enter data copyin(activeTracersCur) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 043ceb8603f4..bcb2a3e01ff6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -1607,7 +1607,7 @@ subroutine ocn_diagnostic_solve_GMvel(layerThickEdgeFlux, normalGMBolusVelocity, real (kind=RKIND), dimension(:), allocatable:: div_huGMBolus #ifdef MPAS_OPENACC - !$acc declare device_resident(div_hu, div_huTransport, div_huGMBolus, div_huMLEBolus) + !$acc declare device_resident(div_huGMBolus) #endif ! Need owned cells for relativeVorticityCell @@ -1619,7 +1619,7 @@ subroutine ocn_diagnostic_solve_GMvel(layerThickEdgeFlux, normalGMBolusVelocity, #ifdef MPAS_OPENACC !$acc parallel loop gang vector & !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & - !$acc areaCell, layerThicknessEdgeFlux, normalGMBolusVelocity, & + !$acc normalGMBolusVelocity, & !$acc vertGMBolusVelocityTop, invAreaCell, minLevelCell) & !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & !$acc dvEdge_temp, k) @@ -1718,7 +1718,7 @@ subroutine ocn_diagnostic_solve_MLEvel(layerThickEdgeFlux, normalMLEVelocity, & #ifdef MPAS_OPENACC !$acc parallel loop gang vector & !$acc present(nEdgesOnCell, edgesOnCell, edgeSignonCell, dcEdge, dvEdge, & - !$acc areaCell, layerThicknessEdgeFlux, normalMLEVelocity, & + !$acc areaCell, normalMLEVelocity, & !$acc vertMLEBolusVelocityTop, invAreaCell, minLevelCell) & !$acc private(invAreaCell1, i, iEdge, edgeSignOnCell_temp, dcEdge_temp, & !$acc dvEdge_temp, k)