diff --git a/components/mpas-ocean/meke_histograms.py b/components/mpas-ocean/meke_histograms.py new file mode 100644 index 000000000000..f81edd6adde4 --- /dev/null +++ b/components/mpas-ocean/meke_histograms.py @@ -0,0 +1,106 @@ +import numpy +import xarray +import matplotlib.pyplot as plt + +comparison=False +if comparison: + master_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-master/ocean/global_ocean/QU240/PHC/performance_test' +#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-off/ocean/global_ocean/QU240/PHC/performance_test' +#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke/ocean/global_ocean/QU240/PHC/performance_test' +meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-rk4/ocean/baroclinic_channel/10km/rpe_test/' +#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-rk4/ocean/baroclinic_channel/10km/default' +#subdir='forward' +subdir='rpe_test_5_nu_200' +if comparison: + master_data = xarray.open_dataset(f'{master_path}/{subdir}/output.nc') +meke_data = xarray.open_dataset(f'{meke_path}/{subdir}/output.nc') + +nt = meke_data.sizes['Time'] + +if comparison: + gmBolusKappa1 = master_data.gmBolusKappa + gmCellsAvg1 = numpy.mean(gmBolusKappa1, axis=1) # nt + gmCellsMin1 = numpy.amin(gmBolusKappa1, axis=1) # nt + gmCellsMax1 = numpy.amax(gmBolusKappa1, axis=1) # nt + +gmBolusKappa = meke_data.gmBolusKappa +gmCellsAvg = numpy.mean(gmBolusKappa, axis=1) # nt +gmCellsMin = numpy.amin(gmBolusKappa, axis=1) # nt +gmCellsMax = numpy.amax(gmBolusKappa, axis=1) # nt + + +mekeSource = meke_data.mekeSourceTopOfEdge.values +mekeSink = meke_data.mekeSink.values +meke = meke_data.meke.values + +mekeCellsAvg = numpy.mean(meke, axis=(1,2)) # nt, nVertLevels +mekeCellsMin = numpy.amin(meke, axis=(1,2)) # nt +mekeCellsMax = numpy.amax(meke, axis=(1,2)) # nt + +mekeSourceEdgesAvg = numpy.mean(mekeSource, axis=(1,2)) # nt, nVertLevels +mekeSourceEdgesMin = numpy.amin(mekeSource, axis=(1,2)) # nt +mekeSourceEdgesMax = numpy.amax(mekeSource, axis=(1,2)) # nt + +mekeSinkCellsAvg = numpy.mean(mekeSink, axis=(1,2)) # nt, nVertLevels +mekeSinkCellsMin = numpy.amin(mekeSink, axis=(1,2)) # nt +mekeSinkCellsMax = numpy.amax(mekeSink, axis=(1,2)) # nt + +tidx = 0 +#tidx = nt-1 + +print('plotting gm') +fig = plt.figure() +if comparison: + plt.hist(gmBolusKappa1[tidx,:],label='master',bins=100) +plt.hist(gmBolusKappa[tidx,:],label='meke',bins=100) +plt.legend() +plt.savefig(f'{meke_path}/gm_histogram.png') +plt.close(fig) + +fig,ax = plt.subplots(2) +print('plotting meke terms') +ax[0].hist(mekeSource[tidx,:].flatten(),label='source',bins=100) +ax[0].hist(mekeSink[tidx,:].flatten(),label='sink',bins=100) +ax[0].set_xlabel('MEKE terms') +ax[0].set_xscale('log') +ax[0].legend() + +print('plotting meke') +ax[1].hist(meke[tidx,:].flatten(),bins=100) +ax[1].set_xlabel('MEKE') +ax[1].set_xscale('log') +ax[1].legend() +fig.savefig(f'{meke_path}/meke_histogram.png') +plt.close(fig) + +fig,ax = plt.subplots(3,1) +ax[0].plot(mekeSourceEdgesMin, '--b') +ax[0].plot(mekeSourceEdgesAvg, '-b', label='Source') +ax[0].plot(mekeSourceEdgesMax, ':b') +ax[0].plot(mekeSinkCellsMin, '--k') +ax[0].plot(mekeSinkCellsAvg, '-k', label='Sink') +ax[0].plot(mekeSinkCellsMax, ':k') +ax[0].set_yscale('log') +ax[0].legend() +ax[1].plot(mekeCellsMin, '--k') +ax[1].plot(mekeCellsAvg, '-k', label='MEKE') +ax[1].plot(mekeCellsMax, ':k') +ax[1].set_yscale('log') +ax[1].legend() +ax[2].plot(gmCellsMin, '--k') +ax[2].plot(gmCellsAvg, '-k', label='GM MEKE') +ax[2].plot(gmCellsMax, ':k') +if comparison: + ax[2].plot(gmCellsMin1, '--b') + ax[2].plot(gmCellsAvg1, '-b', label='GM master') + ax[2].plot(gmCellsMax1, ':b') + ax[2].legend() +fig.savefig(f'{meke_path}/meke_timeseries.png') +plt.close(fig) + +print('mekeSource: min, 10th, 50th, 90th, max') +print(f'{numpy.min(mekeSource[tidx,:])},{numpy.percentile(mekeSource[tidx,:],[10,50,90])},{numpy.max(mekeSource[tidx,:])}') +print('mekeSink: min, 10th, 50th, 90th, max') +print(f'{numpy.min(mekeSink[tidx,:])},{numpy.percentile(mekeSink[tidx,:],[10,50,90])},{numpy.max(mekeSink[tidx,:])}') +print('meke: min, 10th, 50th, 90th, max') +print(f'{numpy.min(meke[tidx,:])},{numpy.percentile(meke[tidx,:],[10,50,90])},{numpy.max(meke[tidx,:])}') diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index fda2d17e651c..eba0a6c1d682 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -348,6 +348,14 @@ description="If true, the standard GM for the tracer advection and mixing is turned on." possible_values=".true. or .false." /> + + + @@ -2129,6 +2138,12 @@ packages="thicknessFilter" /> + + + + + @@ -2853,6 +2875,14 @@ description="gradient Richardson number defined at the edge (horizontally) and top (vertically)" packages="forwardMode;analysisMode" /> + + 0) then @@ -555,6 +579,22 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #endif endif + if (associated(mekeNew)) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k = 1,nVertLevels + mekeNew(k,iCell) = & + mekeCur(k,iCell) + end do + end do + !$omp end do + !$omp end parallel + iCell = 10 + write(log_string,*) 'SE prep: mekeCur, mekeNew', mekeCur(1,iCell), mekeNew(1,iCell) + call mpas_log_write(log_string) + end if + call mpas_timer_stop("se prep") call mpas_pool_get_array(forcingPool, 'seaIcePressure', seaIcePressure) @@ -1754,8 +1794,42 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ end do call mpas_timer_stop("se halo tracers") + if (MEKEActive) then + !if (config_prognostic_meke) then + call mpas_timer_start('se meke tend') + call ocn_tend_meke(tendPool, statePool, & + meshPool, dt, 2) + call mpas_timer_stop('se meke tend') + + ! update halo for tracer tendencies + call mpas_timer_start("se halo tracers") + call mpas_dmpar_field_halo_exch(domain, groupItr%memberName) + call mpas_timer_stop("se halo tracers") + end if + call mpas_timer_start('se loop fini') + call mpas_pool_get_array(tracersPool, 'activeTracers', & + tracersGroupCur, 1) + call mpas_pool_get_array(tracersPool, 'activeTracers', & + tracersGroupNew, 2) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocityCur, 1) + if (MEKEActive) then + call mpas_pool_get_array(statePool, 'meke', & + mekeCur, 1) + call mpas_pool_get_array(statePool, 'meke', & + mekeNew, 2) + call mpas_pool_get_array(tendPool, 'meke', & + mekeTend) + iCell = 10 + write(log_string,*) 'SE after tend: mekeCur, mekeNew', mekeCur(1,iCell), mekeNew(1,iCell) + write(log_string,*) 'SE after tend: mekeTend', mekeTend(1,iCell) + call mpas_log_write(log_string) + endif + call mpas_pool_get_array(tracersTendPool,'activeTracersTend', & + activeTracersTend) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! If iterating, reset variables for next iteration !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1851,6 +1925,39 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #endif end if + if (MEKEActive) then + !if (config_prognostic_meke) then + + !$omp parallel + !$omp do schedule(runtime) private(i, k, temp_h, temp) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + + ! this is h_{n+1} + temp_h = layerThicknessCur(k,iCell) + dt* & + layerThicknessTend(k,iCell) + + ! this is h_{n+1/2} + layerThicknessNew(k,iCell) = 0.5* & + (layerThicknessCur(k,iCell) + temp_h) + + ! This is Phi at n+1 + temp = (mekeCur(k,iCell)* & + layerThicknessCur(k,iCell) + dt* & + mekeTend(k,iCell))/temp_h + + ! This is Phi at n+1/2 + mekeNew(k,iCell) = max(0.0_RKIND, 0.5_RKIND* & + (mekeCur(k,iCell) + temp)) + end do ! vertical + end do ! iCell + !$omp end do + !$omp end parallel + iCell = 10 + write(log_string,*) 'SE update: mekeCur, mekeNew', mekeCur(1,iCell), mekeNew(1,iCell) + call mpas_log_write(log_string) + end if + #ifdef MPAS_OPENACC !$acc parallel loop collapse(2) present(normalVelocityNew) & !$acc present(normalBaroclinicVelocityNew, edgeMask) & @@ -2188,6 +2295,25 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ #endif end if + if (MEKEActive) then + !if (config_prognostic_meke) then + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + mekeNew(k,iCell) = max(0.0_RKIND, & + (mekeCur(k,iCell) * & + layerThicknessCur(k,iCell) + dt* & + mekeTend(k,iCell))/ & + layerThicknessNew(k,iCell)) + end do + end do + !$omp end do + !$omp end parallel + iCell = 10 + write(log_string,*) 'SE update division: mekeCur, mekeNew', mekeCur(1,iCell), mekeNew(1,iCell) + call mpas_log_write(log_string) + end if ! Recompute final u to go on to next step. ! u_{n+1} = normalBarotropicVelocity_{n+1} + diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 81d8235bb707..41a06eb88cf6 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -48,6 +48,7 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_tracer_exponential_decay.o \ mpas_ocn_tracer_ideal_age.o \ mpas_ocn_tracer_TTD.o \ + mpas_ocn_tracer_meke_tend.o \ mpas_ocn_tracer_ecosys.o \ mpas_ocn_tracer_DMS.o \ mpas_ocn_tracer_MacroMolecules.o \ @@ -78,7 +79,7 @@ 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_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_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_meke_tend.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 @@ -96,7 +97,7 @@ 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_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o @@ -140,6 +141,8 @@ mpas_ocn_tracer_hmix_redi.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_high_freq_thickness_hmix_del2.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o +mpas_ocn_tracer_meke_tend.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o + mpas_ocn_tracer_nonlocalflux.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_tracer_short_wave_absorption.o: mpas_ocn_tracer_short_wave_absorption_jerlov.o mpas_ocn_tracer_short_wave_absorption_variable.o mpas_ocn_constants.o mpas_ocn_config.o 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..03b0b8095476 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -54,6 +54,8 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:), pointer :: pressure real (kind=RKIND), dimension(:,:), pointer :: thermExpCoeff real (kind=RKIND), dimension(:,:), pointer :: salineContractCoeff + real (kind=RKIND), dimension(:,:), pointer :: mekeSourceTopOfEdge + real (kind=RKIND), dimension(:,:), pointer :: mekeSink real (kind=RKIND), dimension(:,:), pointer :: BruntVaisalaFreqTop real (kind=RKIND), dimension(:,:), pointer :: tangentialVelocity real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeFlux @@ -142,6 +144,7 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:,:), pointer :: activeTracerHorMixTendency real (kind=RKIND), dimension(:,:), pointer :: temperatureShortWaveTendency real (kind=RKIND), dimension(:), pointer :: RediKappa + real (kind=RKIND), dimension(:), pointer :: LrEdge ! pointers for tendencies used in diagnostic budget computation real (kind=RKIND), dimension(:,:,:), pointer, contiguous :: & @@ -263,8 +266,15 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e gmHorizontalTaper) call mpas_pool_get_array(diagnosticsPool, 'RossbyRadius', & RossbyRadius) + call mpas_pool_get_array(diagnosticsPool, 'LrEdge', & + LrEdge) end if + if ( config_prognostic_meke ) then + call mpas_pool_get_array(diagnosticsPool, 'mekeSourceTopOfEdge', mekeSourceTopOfEdge) + call mpas_pool_get_array(diagnosticsPool, 'mekeSink', mekeSink) + endif + if ( config_compute_active_tracer_budgets ) then call mpas_pool_get_array(diagnosticsPool, & 'activeTracerSurfaceFluxTendency', & @@ -747,6 +757,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e if (config_use_GM.or.config_use_Redi) then !$acc enter data create( & !$acc RediKappa, & + !$acc LrEdge, & !$acc gmBolusKappa, & !$acc gmKappaScaling, & !$acc RediKappaScaling, & @@ -987,6 +998,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ if (config_use_GM.or.config_use_Redi) then !$acc exit data delete( & !$acc RediKappa, & + !$acc LrEdge, & !$acc gmBolusKappa, & !$acc gmKappaScaling, & !$acc RediKappaScaling, & @@ -1181,6 +1193,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ end if if (config_use_GM.or.config_use_Redi) then nullify(RediKappa, & + LrEdge, & gmBolusKappa, & gmKappaScaling, & RediKappaScaling, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index 0e6cdca08cec..6d66754bb428 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -15,10 +15,11 @@ module ocn_gm use ocn_constants use ocn_config + use ocn_mesh use ocn_diagnostics_variables implicit none - private + public save !-------------------------------------------------------------------- @@ -121,6 +122,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND), dimension(:, :), pointer :: & normalVelocity, & + meke, & layerThickness real(kind=RKIND), dimension(:), pointer :: & @@ -135,8 +137,9 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND) :: dcEdgeInv, drhoDx, drhoDT, drhoDS, dTdx, dSdx real(kind=RKIND) :: slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, invAreaCell real(kind=RKIND) :: lt1, lt2, c_min - real(kind=RKIND) :: sigma, Lr, Length, L_rhines, shearEdgeInv - real (kind=RKIND) :: eddyLength, eddyTime + real(kind=RKIND) :: sigma, Length, L_rhines, mekeEdge, shearEdgeInv, Velocity + real(kind=RKIND) :: eddyLength, eddyTime + real(kind=RKIND) :: gmBolusKappaDia real(kind=RKIND), dimension(:), allocatable :: dzTop, dTdzTop, dSdzTop, k33Norm real(kind=RKIND), dimension(:), allocatable :: RossbyRadiusTaper real(kind=RKIND) :: c_Visbeck ! baroclinic wave speed from Visbeck parameterization @@ -156,6 +159,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & real(kind=RKIND), dimension(:, :, :), pointer :: activeTracers integer, pointer :: indexTemperature, indexSalinity integer :: timeLevel + character(len=100) :: log_string if (present(timeLevelIn)) then timeLevel = timeLevelIn @@ -169,6 +173,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) + call mpas_pool_get_array(statePool, 'meke', meke, timeLevel) call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, timeLevel) call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) @@ -652,35 +657,68 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & if (local_config_GM_compute_EdenGreatbatch) then call mpas_pool_get_array(meshPool, 'fEdge', fEdge) - !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, sigma, RiTopOfEdge, & - !$omp L_rhines, Lr, Length, BruntVaisalaFreqTopEdge, shearEdgeInv) - do iEdge=1,nEdges - cell1 = cellsOnEdge(1, iEdge) - cell2 = cellsOnEdge(2, iEdge) - Lr = min(cGMphaseSpeed(iEdge) / (1.0E-15_RKIND + abs(fEdge(iEdge))), & - sqrt(cGMphaseSpeed(iEdge) / (2.0_RKIND*betaEdge(iEdge)))) - - do k=minLevelEdgeBot(iEdge)+1,maxLevelEdgeTop(iEdge) - BruntVaisalaFreqTopEdge = 0.5_RKIND*(max(BruntVaisalaFreqTop(k,cell1),0.0_RKIND) + & - max(BruntVaisalaFreqTop(k,cell2),0.0_RKIND)) - shearEdgeInv = (0.5_RKIND*(layerThickEdgeMean(k-1,iEdge) + layerThickEdgeMean(k,iEdge)))**2.0 & - / ( (normalVelocity(k-1,iEdge) - normalVelocity(k,iEdge))**2.0 & - + (tangentialVelocity(k-1,iEdge) - tangentialVelocity(k,iEdge))**2.0 & - + 1.0e-18_RKIND ) - RiTopOfEdge = BruntVaisalaFreqTopEdge * shearEdgeInv - sigma = max(abs(fEdge(iEdge)),sqrt(2.0_RKIND*betaEdge(iEdge)* & - cGMphaseSpeed(iEdge))) / sqrt(RiTopOfEdge + config_GM_EG_riMin) - L_rhines = sigma / (1.0E-18_RKIND + betaEdge(iEdge)) - Length = min(config_GM_EG_Rossby_factor*Lr,config_GM_EG_Rhines_factor*L_rhines) - !For compatibility with other schemes, we normalize by the 2-D kappa field, which will - !be removed below in the computation of the stream function - gmKappaScaling(k,iEdge) = max(config_GM_spatially_variable_min_kappa, min(config_GM_EG_kappa_factor* & - Length**2.0*sigma, config_GM_spatially_variable_max_kappa)) / gmBolusKappa(iEdge) + if (config_prognostic_meke .or. config_use_prognostic_meke) then + !$omp parallel + !$omp do schedule(runtime) + do iEdge = 1,nEdges + if (sphereRadius == 0.0_RKIND) then + LrEdge(iEdge) = cGMphaseSpeed(iEdge) / (1.0E-15_RKIND + abs(fEdge(iEdge))) + else + LrEdge(iEdge) = min(cGMphaseSpeed(iEdge) / (1.0E-15_RKIND + abs(fEdge(iEdge))), & + sqrt(cGMphaseSpeed(iEdge) / (2.0_RKIND*betaEdge(iEdge)))) + endif enddo - enddo - !$omp end do - !$omp end parallel + !$omp end do + !$omp end parallel + endif + if (config_use_prognostic_meke) then + !$omp parallel + !$omp do schedule(runtime) private(cell1, cell2, & + !$omp Length, mekeEdge, Velocity) + do iEdge = 1,nEdges + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + call mpas_log_write('Assigning gmKappaScaling based on meke') + do k = minLevelEdgeBot(iEdge)+1,maxLevelEdgeTop(iEdge) + mekeEdge = 0.5_RKIND*(meke(k,cell1) + meke(k,cell2)) + Velocity = mekeEdge**0.5_RKIND + Length = min((Velocity/betaEdge(iEdge))**0.5_RKIND, LrEdge(iEdge)) + gmKappaScaling(k,iEdge) = max(config_GM_spatially_variable_min_kappa, & + min(config_GM_spatially_variable_max_kappa, & + Length*Velocity) & + ) / gmBolusKappa(iEdge) + enddo + enddo + !$omp end do + !$omp end parallel + else + !$omp parallel + !$omp do schedule(runtime) private(cell1, cell2, sigma, RiTopOfEdge, & + !$omp L_rhines, Length, BruntVaisalaFreqTopEdge, shearEdgeInv) + do iEdge=1,nEdges + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + do k=minLevelEdgeBot(iEdge)+1,maxLevelEdgeTop(iEdge) + BruntVaisalaFreqTopEdge = 0.5_RKIND*(max(BruntVaisalaFreqTop(k,cell1),0.0_RKIND) + & + max(BruntVaisalaFreqTop(k,cell2),0.0_RKIND)) + shearEdgeInv = (0.5_RKIND*(layerThickEdgeMean(k-1,iEdge) + layerThickEdgeMean(k,iEdge)))**2.0 & + / ( (normalVelocity(k-1,iEdge) - normalVelocity(k,iEdge))**2.0 & + + (tangentialVelocity(k-1,iEdge) - tangentialVelocity(k,iEdge))**2.0 & + + 1.0e-18_RKIND ) + RiTopOfEdge = BruntVaisalaFreqTopEdge * shearEdgeInv + sigma = max(abs(fEdge(iEdge)),sqrt(2.0_RKIND*betaEdge(iEdge)* & + cGMphaseSpeed(iEdge))) / sqrt(RiTopOfEdge + config_GM_EG_riMin) + L_rhines = sigma / (1.0E-18_RKIND + betaEdge(iEdge)) + Length = min(config_GM_EG_Rossby_factor*LrEdge(iEdge),config_GM_EG_Rhines_factor*L_rhines) + !For compatibility with other schemes, we normalize by the 2-D kappa field, which will + !be removed below in the computation of the stream function + gmKappaScaling(k,iEdge) = max(config_GM_spatially_variable_min_kappa, min(config_GM_EG_kappa_factor* & + Length**2.0*sigma, config_GM_spatially_variable_max_kappa)) / gmBolusKappa(iEdge) + enddo + enddo + !$omp end do + !$omp end parallel + endif endif nEdges = nEdgesArray(3) @@ -864,6 +902,35 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & deallocate (tridiagB) deallocate (tridiagC) + if (config_prognostic_meke) then + !$omp parallel + !$omp do schedule(runtime) private(k, cell1, cell2, BruntVaisalaFreqTopEdge) + do iEdge = 1, nEdges + + mekeSourceTopOfEdge(:, iEdge) = 0.0_RKIND + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + BruntVaisalaFreqTopEdge = 0.5_RKIND*(BruntVaisalaFreqTop(k, cell1) + BruntVaisalaFreqTop(k, cell2)) + BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) + if (BruntVaisalaFreqTopEdge < 1.0E-10_RKIND) then + gmBolusKappaDia = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge) + else + gmBolusKappaDia = 0.0_RKIND + endif + mekeSourceTopOfEdge(k, iEdge) = ( gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge) - & + gmBolusKappaDia ) * & + (gradDensityConstZTopOfEdge(k,iEdge))**2 / & + (BruntVaisalaFreqTopEdge + 1.0E-10_RKIND)! & + !- gmBolusKappaDia*BruntVaisalaFreqTopEdge + !> I can't figure out why this term would be negative + end do + end do + !$omp end do + !$omp end parallel + end if + end if !end config_use_GM deallocate(gradDensityEdge, & dDensityDzTopOfCell, & @@ -968,6 +1035,8 @@ subroutine ocn_GM_init(domain, err)!{{{ real(kind=RKIND), pointer :: sphere_radius real(kind=RKIND), dimension(:), pointer :: latEdge, fEdge integer, dimension(:, :), pointer :: cellsOnEdge + real(kind=RKIND), dimension(:), allocatable :: Lr + character(len=100) :: log_string err = 0 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F index 65b839b1005b..b26407fe4b77 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_init_routines.F @@ -105,6 +105,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ integer, dimension(:,:), pointer :: boundaryCellTmp real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:,:), pointer :: normalVelocity + real (kind=RKIND), dimension(:,:), pointer :: meke real (kind=RKIND), dimension(:,:,:), pointer :: tracersGroup @@ -121,6 +122,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers type (mpas_pool_iterator_type) :: groupItr + character (len=100) :: log_string call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) @@ -145,6 +147,7 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure) call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1) + call mpas_pool_get_array(statePool, 'meke', meke, 1) if (.not. config_do_restart) then do iCell=1,nCells @@ -234,6 +237,13 @@ subroutine ocn_init_routines_block(block, dt, err)!{{{ end if end do + ! Temporarily set equal to temperature + if (config_prognostic_meke) then + do iCell=1,nCells + meke(:, iCell) = 0.0_RKIND + end do + end if + ! ------------------------------------------------------------------ ! Accumulating various parametrizations of the transport velocity ! ------------------------------------------------------------------ diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index e081df8ae097..79ce2351b6eb 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -51,6 +51,7 @@ module ocn_tendency use ocn_tracer_surface_flux_to_tend use ocn_tracer_ecosys use ocn_tracer_DMS + use ocn_tracer_meke_tend use ocn_tracer_MacroMolecules use ocn_tracer_CFC @@ -88,6 +89,7 @@ module ocn_tendency ocn_tend_vel, & ocn_tend_tracer, & ocn_tend_freq_filtered_thickness, & + ocn_tend_meke, & ocn_tendency_init !-------------------------------------------------------------------- @@ -1306,6 +1308,199 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & end subroutine ocn_tend_tracer!}}} +!*********************************************************************** +! +! routine ocn_tend_meke +! +!> \brief Computes mesoscale eddy kinetic energy tendency +!> \author Carolyn Begeman, Luke Van Roekel +!> \date June 2022 +!> \details +!> This routine computes meke tendency for the ocean +! +!----------------------------------------------------------------------- + + subroutine ocn_tend_meke(tendPool, statePool, & + meshPool, dt, timeLevelIn )!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + statePool, &!< [in] ocean state variables + meshPool !< [in] mesh information + + real (kind=RKIND), intent(in) :: & + dt !< [in] time step (seconds) + + integer, intent(in), optional :: & + timeLevelIn !< [in] time index to use for state vars + + ! these variables input from shared diagnostics module: + ! normalTransportVelocity - transport velocity across edge + ! layerThickEdgeFlux - flux-related layer thickness on edge + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(inout) :: & + tendPool !< [inout] Tendency terms for all variables + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + real (kind=RKIND), dimension(:,:), pointer, contiguous :: & + meke, &! mesoscale EKE + mekeTend ! mesoscale EKE tendency + + real (kind=RKIND), dimension(:,:,:), allocatable :: & + meke3d, &! mesoscale EKE + mekeTend3d ! mesoscale EKE tendency + + integer :: & + iCell, iEdge, k, &! loop counters for cell,edge,vert,tracer + err, &! internal error flag + timeLevel ! time level to use for state variables + + real (kind=RKIND), dimension(:,:), pointer, contiguous :: & + layerThickness ! layer thickness + + ! Scratch Arrays + real (kind=RKIND), dimension(:,:), allocatable :: & + normalThicknessFlux ! Flux of thickness through edge (m^2/s) + + character(len=100) :: log_string + ! these variables input from shared diagnostics module: + ! normalTransportVelocity - transport velocity across edge + ! layerThickEdgeFlux - flux-related layer thickness on edge + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + call mpas_log_write("ocn_tend_meke") + call mpas_timer_start("ocn_tend_meke") + + ! set time level to default 1 or override with input arg + if (present(timeLevelIn)) then + timeLevel = timeLevelIn + else + timeLevel = 1 + end if + + ! retrieve data from pools + call mpas_pool_get_array(statePool, 'meke', & + meke, timeLevel) + call mpas_pool_get_array(tendPool, 'meke', & + mekeTend) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThickness, timeLevel) + iCell = 10 + write(log_string,*) 'meke, mekeTend =',meke(1,iCell),mekeTend(1,iCell) + call mpas_log_write(log_string) + + allocate(meke3d (1, nVertLevels, nCellsAll+1), & + mekeTend3d(1, nVertLevels, nCellsAll+1)) + + ! initialize meke tendencies to zero. + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k=1, nVertLevels + mekeTend3d(1, k, iCell) = 0.0_RKIND + end do + end do + !$omp end do + !$omp end parallel + + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k=1, nVertLevels + meke3d(1, k, iCell) = meke(k, iCell) + end do + end do + !$omp end do + !$omp end parallel + + ! allocate and transfer data not specific to tracer groups + allocate(normalThicknessFlux(nVertLevels, nEdgesAll+1)) + + !$acc enter data create(normalThicknessFlux) & + !$acc copyin(layerThickness) + !TEMPORARY - once diagnotics module ported, these should already + ! be on the device + !$acc update device (normalTransportVelocity, layerThickEdgeFlux, & + !$acc vertAleTransportTop) + + ! compute transport velocity to use for all tracers +#ifdef MPAS_OPENACC + !$acc parallel loop collapse(2) & + !$acc present(normalThicknessFlux, normalTransportVelocity, & + !$acc layerThickEdgeFlux) +#else + !$omp parallel + !$omp do schedule(runtime) private(k) +#endif + do iEdge = 1,nEdgesAll + do k = 1,nVertLevels + normalThicknessFlux(k,iEdge) = & + normalTransportVelocity(k,iEdge)* & + layerThickEdgeFlux(k, iEdge) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + ! Accumulate tendencies for all tracers + ! First is the horizontal advection term + ! -div( layerThickness \phi u) + ! Tracer budget is computed and stored within the tracer adv + ! routine + call mpas_log_write('Advect MEKE') + call ocn_tracer_advection_tend(meke3d, & + normalThicknessFlux, vertAleTransportTop, & + layerThickness, dt, mekeTend3d, .false.) + call mpas_log_write('... completed') + + ! mekeTend3d has an additional dimension, which can hold + ! multiple tracers but doesn't in the case of MEKE + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsAll + do k=1, nVertLevels + mekeTend(k, iCell) = mekeTend3d(1, k, iCell) + end do + end do + !$omp end do + !$omp end parallel + iCell = 10 + write(log_string,*) 'After adv mekeTend =',mekeTend(1,iCell) + call mpas_log_write(log_string) + + call mpas_log_write('MEKE source') + call ocn_tracer_meke_source_tend(layerThickness, dt, mekeTend, err) + call mpas_log_write('... completed') + iCell = 10 + write(log_string,*) 'After source mekeTend =',mekeTend(1,iCell) + call mpas_log_write(log_string) + + call mpas_log_write('MEKE sink') + call ocn_tracer_meke_sink_tend(meke, layerThickness, dt, mekeTend, err) + call mpas_log_write('... completed') + iCell = 10 + write(log_string,*) 'After sink meke, mekeTend =',meke(1,iCell),mekeTend(1,iCell) + call mpas_log_write(log_string) + + call mpas_timer_stop("ocn_tend_meke") + + !-------------------------------------------------------------------- + + end subroutine ocn_tend_meke!}}} + !*********************************************************************** ! ! routine ocn_tend_freq_filtered_thickness diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F index b2fe794421af..b88a27e0257e 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F @@ -141,6 +141,7 @@ subroutine ocn_tracer_advection_init(err)!{{{ err = 0 ! initialize error code to success + call mpas_log_write('Start tracer adv init') ! set some basic flags for options tracerAdvOff = config_disable_tr_adv if (config_flux_limiter == 'monotonic') then @@ -153,10 +154,15 @@ subroutine ocn_tracer_advection_init(err)!{{{ endif ! set all other options from submodule initialization routines + call mpas_log_write('Start tracer shared init') call ocn_tracer_advection_shared_init(err1) + call mpas_log_write('End shared init') call ocn_tracer_advection_std_init (err2) + call mpas_log_write('End std init') call ocn_tracer_advection_mono_init (err3) + call mpas_log_write('End mono init') call ocn_tracer_advection_vert_init (err4) + call mpas_log_write('End vert init') ! if an error is returned from init routines, write an error ! message and return a non-zero error code diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_meke_tend.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_meke_tend.F new file mode 100644 index 000000000000..a712536b7b40 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_meke_tend.F @@ -0,0 +1,317 @@ +! 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 +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_tracer_meke +! +!> \brief MPAS ocean Mesoscale EKE source and sink terms +!> \author Carolyn Begeman +!> \date 07/11/2022 +!> \version +!> \details +!> This module contains the routines for computing MEKE tendencies +!> due to source and sink terms +! +!----------------------------------------------------------------------- + +module ocn_tracer_meke_tend + + use mpas_timer + use mpas_derived_types + use mpas_pool_routines + use ocn_constants + use ocn_config + use ocn_diagnostics_variables + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_tracer_meke_source_tend, & + ocn_tracer_meke_sink_tend, & + ocn_tracer_meke_tend_init + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_tracer_meke_source_tend +! +!> \brief Computes MEKE source terms +!> \author Carolyn Begeman +!> \date 07/11/2022 +!> \details +!> This routine computes the tendency for MEKE due to baroclinic +!> instability source +! +!----------------------------------------------------------------------- + + subroutine ocn_tracer_meke_source_tend(layerThickness, dt, tend, err)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), intent(in) :: dt ! time step + real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: MEKE tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, iEdge, k, nCells, nEdges + real (kind=RKIND), dimension(:,:), allocatable :: & + mekeSource, & + mekeSourceEdge + + call mpas_log_write('MEKE source start') + err = 0 + + if (.not. config_prognostic_meke) return + + call mpas_timer_start('MEKE source') + + nCells = nCellsOwned + nEdges = nEdgesOwned + + allocate(mekeSource(nVertLevels, nCells)) + allocate(mekeSourceEdge(nVertLevels+1, nEdges)) + + ! mekeSourceTopOfEdge comes from diagnostics_variables + ! first, interpolate mekeSourceTopOfEdge to mekeSourceEdge + !$omp parallel + !$omp do schedule(runtime) private(k) + do iEdge = 1, nEdges + mekeSourceEdge(:, iEdge) = 0.0_RKIND + do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge) + mekeSourceEdge(k, iEdge) = 0.5_RKIND * & + (mekeSourceTopOfEdge(k, iEdge) + & + mekeSourceTopOfEdge(k+1, iEdge)) + end do + end do + !$omp end do + !$omp end parallel + + ! second, interpolate mekeSourceEdge to cell centers + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + mekeSource(:, iCell) = 0.0_RKIND + do iEdge = 1, nEdgesOnCell(iCell) + do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + mekeSource(k, iCell) = mekeSource(k, iCell) + & + mekeSourceEdge(k, iEdge) * edgeAreaFractionOfCell(iEdge, iCell) + end do + end do + end do + !$omp end do + !$omp end parallel + + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + do k = minLevelCell(iCell), maxLevelCell(iCell) + + tend(k, iCell) = tend(k, iCell) & + + dt * mekeSource(k, iCell) * layerThickness(k, iCell) + end do + end do + !$omp end do + !$omp end parallel + + call mpas_timer_stop('MEKE source') + + !-------------------------------------------------------------------- + + end subroutine ocn_tracer_meke_source_tend!}}} + +!*********************************************************************** +! +! routine ocn_tracer_meke_sink_tend +! +!> \brief Computes MEKE sink terms +!> \author Carolyn Begeman +!> \date 07/11/2022 +!> \details +!> This routine computes the tendency for MEKE due to dissipation +! +!----------------------------------------------------------------------- + + subroutine ocn_tracer_meke_sink_tend(meke, layerThickness, dt, tend, err)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: meke ! MEKE + real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness + real (kind=RKIND), intent(in) :: dt ! time step + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: MEKE tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: cell1, cell2, iCell, iEdge, j, k, nCells, nEdges + real(kind=RKIND), parameter :: c_diss = 0.1_RKIND + real(kind=RKIND) :: Lr, Velocity, Length, betaCell + character(len=100) :: log_string + + call mpas_log_write('MEKE sink start') + iCell = 10 + write(log_string,*) 'meke =',meke(1,iCell) + call mpas_log_write(log_string) + err = 0 + + if (.not. config_prognostic_meke) return + + call mpas_timer_start('MEKE sink') + + nCells = nCellsOwned + nEdges = nEdgesOwned + + !$omp parallel + !$omp do schedule(runtime) private(iEdge, j, k, Lr, betaCell, Velocity, Length) + do iCell=1,nCells + Lr = 0.0_RKIND + do j = 1,nEdgesOnCell(iCell) + iEdge = edgesOnCell(j,iCell) + Lr = Lr + LrEdge(iEdge)*edgeAreaFractionOfCell(j, iCell) + enddo + if (sphereRadius == 0.0_RKIND) then + betaCell = 1e-20_RKIND + else + betaCell = 2.0_RKIND*omega*cos(latCell(iCell)) / sphereRadius + endif + do k = minLevelCell(iCell),maxLevelCell(iCell) + Velocity = meke(k,iCell)**0.5_RKIND + if (sphereRadius == 0.0_RKIND) then + Length = Lr + else + Length = min((Velocity/betaCell)**0.5_RKIND, Lr) + endif + Length = max(1e-10_RKIND, Length) + mekeSink(k,iCell) = c_diss * & + (meke(k,iCell)/layerThickness(k,iCell))**1.5_RKIND / & + Length + if (mekeSink(k,iCell) > 1.0_RKIND .or. (iCell == 10 .and. k==1)) then + write(log_string,*) 'Lr', Lr + call mpas_log_write(log_string) + write(log_string,*) 'omega, sphereRadius', omega, sphereRadius + call mpas_log_write(log_string) + write(log_string,*) 'betaCell, Velocity', betaCell, Velocity + call mpas_log_write(log_string) + write(log_string,*) 'Length, meke', Length, meke(k,iCell) + call mpas_log_write(log_string) + write(log_string,*) 'mekeSink', mekeSink(k,iCell) + call mpas_log_write(log_string) + endif + tend(k, iCell) = tend(k, iCell) & + - dt * mekeSink(k, iCell) * layerThickness(k, iCell) + enddo + enddo + !$omp end do + !$omp end parallel + iCell = 10 + write(log_string,*) 'mekeTend =',tend(1,iCell) + call mpas_log_write(log_string) + call mpas_timer_stop('MEKE sink') + call mpas_log_write('MEKE sink end') + + end subroutine ocn_tracer_meke_sink_tend + +!*********************************************************************** +! +! routine ocn_tracer_meke_tend_init +! +!> \brief Initializes ocean MEKE tendency computation +!> \author Carolyn Begeman +!> \date 07/11/2022 +!> \version +!> \details +!> This routine initializes quantities related to MEKE tendency +! +!----------------------------------------------------------------------- + + subroutine ocn_tracer_meke_tend_init(err)!{{{ + + !-------------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + err = 0 + + end subroutine ocn_tracer_meke_tend_init!}}} + +!*********************************************************************** + +end module ocn_tracer_meke_tend + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker