From dea07e06ed3848aa37f5fe890a605d9c9d0103a5 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 20 Jun 2019 07:01:06 -0700 Subject: [PATCH 1/3] Adds options for 3d varying GM bolus and 2d varying phase speed In the GM bolus calculation, there is a specified diffusivity and phase speed (see Ferrari et al. 2010). Currently MPAS-O assumes these values are fixed in time and space. This PR allows separate flags to enable 2D varying (+time) phase speed (based on the first baroclinic mode phase speed) and a separate flag for 3D (+time) varying diffusivity. --- Registry.xml | 25 ++++++++++ shared/mpas_ocn_gm.F | 107 ++++++++++++++++++++++++++++++++++++------- 2 files changed, 116 insertions(+), 16 deletions(-) diff --git a/Registry.xml b/Registry.xml index 5e5fca2daa13..b4a2c9098894 100644 --- a/Registry.xml +++ b/Registry.xml @@ -375,6 +375,22 @@ description="If true, the standard GM for the tracer advection and mixing is turned on. This includes the redi portion which is captured in hmix." possible_values=".true. or .false." /> + + + + + + 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/countN2)*bottomAv / pii) + + enddo + !$omp end do + + endif + + !$omp do schedule(runtime) + do iEdge=1,nEdges + kappaGM3D(:,iEdge) = gmBolusKappa(iEdge) + enddo + !$omp end do + + if (config_gm_kappa_lat_depth_variable) then + + !$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, maxN) + do iEdge = 1,nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + maxN = -1.0_RKIND + do k=2,maxLevelEdgeTop(iEdge) + BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) + BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) + + maxN = max(maxN,BruntVaisalaFreqTopEdge) + enddo + + do k=2,maxLevelEdgeTop(iEdge) + BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) + BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) + + kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, & + BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND)) + enddo + enddo + !$omp end do + + endif nEdges = nEdgesArray( 3 ) @@ -533,34 +603,34 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool) k = 2 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) - tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1,iEdge) & + tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1,iEdge) & * layerThicknessEdge(k,iEdge)) - BruntVaisalaFreqTopEdge - tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) & + tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) - rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Second to next to the last rows do k = 3, maxLevelEdgeTop(iEdge)-1 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) - tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1, iEdge) & + tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) - tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) & + tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge) ) - BruntVaisalaFreqTopEdge - tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) & + tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) - rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) end do ! Last row k = maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) - tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1,iEdge) & + tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) - tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) & + tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge)) - BruntVaisalaFreqTopEdge - rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) + rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Total number of rows N = maxLevelEdgeTop(iEdge) - 1 @@ -747,7 +817,11 @@ subroutine ocn_gm_init(domain, err)!{{{ call mpas_pool_get_config(ocnConfigs, 'config_use_Redi_bottom_layer_tapering',config_use_Redi_bottom_layer_tapering) call mpas_pool_get_config(ocnConfigs, 'config_Redi_surface_layer_tapering_extent',config_Redi_surface_layer_tapering_extent) call mpas_pool_get_config(ocnConfigs, 'config_Redi_bottom_layer_tapering_depth',config_Redi_bottom_layer_tapering_depth) - + call mpas_pool_get_config(ocnConfigs, 'config_gm_lat_variable_c2',config_gm_lat_variable_c2) + call mpas_pool_get_config(ocnConfigs, 'config_gm_kappa_lat_depth_variable', config_gm_kappa_lat_depth_variable) + call mpas_pool_get_config(ocnConfigs, 'config_gm_kappa_lat_variable', config_gm_kappa_lat_variable) + call mpas_pool_get_config(ocnConfigs, 'config_gm_min_stratification_ratio', config_gm_min_stratification_ratio) + call mpas_pool_get_config(ocnConfigs, 'config_gm_min_phase_speed', config_gm_min_phase_speed) block => domain % blocklist do while (associated(block)) @@ -788,6 +862,7 @@ subroutine ocn_gm_init(domain, err)!{{{ call mpas_dmpar_finalize(domain % dminfo) end if + block => block % next end do end subroutine ocn_gm_init!}}} From dab5e06e4e4cffd0b9e1d03e9b682a9adfb9af7d Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Mon, 12 Aug 2019 19:07:10 -0700 Subject: [PATCH 2/3] Adds 2d variable GM and fixes phase speed bug --- Registry.xml | 8 +++++++ shared/mpas_ocn_gm.F | 51 +++++++++++++++++++++++++++++++++++--------- 2 files changed, 49 insertions(+), 10 deletions(-) diff --git a/Registry.xml b/Registry.xml index b4a2c9098894..98df4df498a4 100644 --- a/Registry.xml +++ b/Registry.xml @@ -383,6 +383,10 @@ description="If true, spatitally and deph varying Bolus Kappa is used, following Danabasoglu et al (2007)" possible_values=".true. or .false." /> + + 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/countN2)*bottomAv / pii) + if(countN2 > 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/ltSum)*ltSum / 3.141592_RKIND) enddo !$omp end do + else + cGMphaseSpeed(:) = config_gravWaveSpeed_trunc endif !$omp do schedule(runtime) @@ -574,6 +577,7 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) maxN = max(maxN,BruntVaisalaFreqTopEdge) + enddo do k=2,maxLevelEdgeTop(iEdge) @@ -582,12 +586,39 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool) kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, & BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND)) + if(kappaGM3D(k,iEdge) <= 0) then + print *, 'val = ',kappaGM3D(k,iEdge),k,iEdge + stop + endif enddo enddo !$omp end do + if (config_gm_kappa_lat_variable) then + !$omp do schedule(runtime) private(cell1, cell2, k, kappaSum, ltSum) + do iEdge = 1, nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + kappaSum = 0.0_RKIND + ltSum = 0.0_RKIND + do k=1,maxLevelEdgeTop(iEdge) + kappaSum = kappaSum + layerThicknessEdge(k,iEdge)*kappaGM3D(k,iEdge) + ltSum = ltSum + layerThicknessEdge(k,iEdge) + enddo + + kappaGM2D(iEdge) = kappaSum / ltSum + !over write the kappa 3d variable + kappaGM3D(:,iEdge) = kappaGM2D(iEdge) + + enddo + !$omp end do + endif + endif + !next do the 2d bit, just average it for now + nEdges = nEdgesArray( 3 ) !$omp do schedule(runtime) From 04874bf7d1d64974983de563afe0a1ea9b9cb829 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Fri, 16 Aug 2019 12:50:09 -0500 Subject: [PATCH 3/3] Removes 2D GM code and debug statements Removes the vertical averaging of the GM kappa bolus diffusivity. Testing shows this functionality is not beneficial to simulations --- Registry.xml | 10 ---------- shared/mpas_ocn_gm.F | 33 +-------------------------------- 2 files changed, 1 insertion(+), 42 deletions(-) diff --git a/Registry.xml b/Registry.xml index 98df4df498a4..2c7f36f454f8 100644 --- a/Registry.xml +++ b/Registry.xml @@ -383,10 +383,6 @@ description="If true, spatitally and deph varying Bolus Kappa is used, following Danabasoglu et al (2007)" possible_values=".true. or .false." /> - - -