Skip to content

Commit

Permalink
Merge PR MPAS-Dev#280 'ocean/redi_add_triads' into ocean/develop
Browse files Browse the repository at this point in the history 
Currently, the Redi along-isopycnal mixing formulation used simple
spatial averaging to compute the isopycnal slope. This PR introduces
triad averaging, which uses a very specific form of an expansion of the
slope using three cells, the thermal expansion and saline contraction
coefficients.
  • Loading branch information
@mark-petersen
mark-petersen committed Jan 27, 2020
2 parents 9558bc3 + 7b60bf6 commit 3e4d773
Show file tree
Hide file tree
Showing 56 changed files with 2,194 additions and 750 deletions.
184 changes: 64 additions & 120 deletions src/core_ocean/Registry.xml
View file

Large diffs are not rendered by default.

18 changes: 5 additions & 13 deletions src/core_ocean/driver/mpas_ocn_core_interface.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,8 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
logical, pointer :: config_use_frazil_ice_formation
logical, pointer :: config_use_tidal_forcing
logical, pointer :: config_use_tidal_potential_forcing
logical, pointer :: config_use_standardGM
logical, pointer :: config_use_GM
logical, pointer :: config_use_Redi
logical, pointer :: config_use_time_varying_atmospheric_forcing
logical, pointer :: config_use_time_varying_land_ice_forcing

Expand Down Expand Up @@ -296,22 +297,13 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
tidalPotentialForcingPKGActive = .true.
end if

!
! test for form of pressure gradient computation
!
! TDR: need to add PKG
call mpas_pool_get_package(packagePool, 'inSituEOSActive', inSituEOSActive)
call mpas_pool_get_config(configPool, 'config_pressure_gradient_type', config_pressure_gradient_type)
if (config_pressure_gradient_type.eq.'Jacobian_from_TS') then
inSituEOSActive = .true.
end if

!
! test for use of gm
!
call mpas_pool_get_package(packagePool, 'gmActive', gmActive)
call mpas_pool_get_config(configPool, 'config_use_standardGM', config_use_standardGM)
if (config_use_standardGM) then
call mpas_pool_get_config(configPool, 'config_use_GM', config_use_GM)
call mpas_pool_get_config(configPool, 'config_use_Redi', config_use_Redi)
if (config_use_GM.or.config_use_Redi) then
gmActive = .true.
end if

Expand Down
2 changes: 1 addition & 1 deletion src/core_ocean/mode_forward/mpas_ocn_forward_mode.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{

call ocn_tracer_hmix_init(err_tmp)
ierr = ior(ierr, err_tmp)
call ocn_tracer_hmix_redi_init(err_tmp)
call ocn_tracer_hmix_redi_init(domain,err_tmp)
ierr = ior(ierr, err_tmp)
call ocn_tracer_surface_flux_init(err_tmp)
ierr = ior(ierr, err_tmp)
Expand Down
28 changes: 15 additions & 13 deletions src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
! Config options
logical, pointer :: config_prescribe_velocity, config_prescribe_thickness
logical, pointer :: config_filter_btr_mode, config_use_freq_filtered_thickness
logical, pointer :: config_use_standardGM
logical, pointer :: config_use_GM
logical, pointer :: config_use_cvmix_kpp
logical, pointer :: config_use_tracerGroup
logical, pointer :: config_disable_thick_all_tend
Expand Down Expand Up @@ -204,7 +204,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call mpas_pool_get_config(domain % configs, 'config_prescribe_velocity', config_prescribe_velocity)
call mpas_pool_get_config(domain % configs, 'config_prescribe_thickness', config_prescribe_thickness)
call mpas_pool_get_config(domain % configs, 'config_use_freq_filtered_thickness', config_use_freq_filtered_thickness)
call mpas_pool_get_config(domain % configs, 'config_use_standardGM', config_use_standardGM)
call mpas_pool_get_config(domain % configs, 'config_use_GM', config_use_GM)
call mpas_pool_get_config(domain % configs, 'config_use_cvmix_kpp', config_use_cvmix_kpp)
call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode)
call mpas_pool_get_config(domain % configs, 'config_disable_vel_all_tend', config_disable_vel_all_tend)
Expand Down Expand Up @@ -726,12 +726,13 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
!$omp end do
! Compute normalGMBolusVelocity and the tracer transport velocity
if (config_use_standardGM) then
call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
if (config_use_GM) then
call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &
meshPool, scratchPool, timeLevelIn=2)
end if
call mpas_threading_barrier()
if (config_use_standardGM) then
if (config_use_GM) then
!$omp do schedule(runtime)
do iEdge = 1, nEdges
normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + normalGMBolusVelocity(:, iEdge)
Expand Down Expand Up @@ -765,7 +766,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
call ocn_time_average_coupled_accumulate(diagnosticsPool, statePool, forcingPool, 2)
if (config_use_standardGM) then
if (config_use_GM) then
call ocn_reconstruct_gm_vectors(diagnosticsPool, meshPool)
end if
call mpas_threading_barrier()
Expand Down Expand Up @@ -1095,7 +1096,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
real (kind=RKIND), intent(in) :: rkSubstepWeight
integer, intent(out) :: err
logical, pointer :: config_prescribe_velocity, config_prescribe_thickness, config_use_standardGM
logical, pointer :: config_prescribe_velocity, config_prescribe_thickness, config_use_GM
integer, pointer :: nCells, nEdges
integer :: iCell, iEdge, k
Expand Down Expand Up @@ -1123,7 +1124,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
call mpas_pool_get_config(block % configs, 'config_prescribe_velocity', config_prescribe_velocity)
call mpas_pool_get_config(block % configs, 'config_prescribe_thickness', config_prescribe_thickness)
call mpas_pool_get_config(block % configs, 'config_use_standardGM', config_use_standardGM)
call mpas_pool_get_config(block % configs, 'config_use_GM', config_use_GM)
call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells)
call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
Expand Down Expand Up @@ -1249,12 +1250,13 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
!$omp end do
! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on
if (config_use_standardGM) then
call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
if (config_use_GM) then
call ocn_gm_compute_Bolus_velocity(provisStatePool, diagnosticsPool, &
meshPool, scratchPool, timeLevelIn=1)
end if
call mpas_threading_barrier()
if (config_use_standardGM) then
if (config_use_GM) then
!$omp do schedule(runtime)
do iEdge = 1, nEdges
normalTransportVelocity(:, iEdge) = normalTransportVelocity(:, iEdge) + normalGMBolusVelocity(:,iEdge)
Expand Down Expand Up @@ -1406,11 +1408,11 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{
character (len=StrKIND) :: modifiedGroupName
character (len=StrKIND) :: configName
logical, pointer :: config_use_standardGM
logical, pointer :: config_use_GM
err = 0
call mpas_pool_get_config(block % configs, 'config_use_standardGM', config_use_standardGM)
call mpas_pool_get_config(block % configs, 'config_use_GM', config_use_GM)
call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells)
call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
Expand Down
85 changes: 73 additions & 12 deletions src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
real (kind=RKIND), dimension(:), allocatable:: uTemp
real (kind=RKIND), dimension(:), pointer :: btrvel_temp
type (field1DReal), pointer :: btrvel_tempField
logical :: activeTracersOnly ! if true only compute tendencies for active tracers
logical :: activeTracersOnly ! if true only compute tendencies for active tracers
integer :: tsIter
integer :: edgeHaloComputeCounter, cellHaloComputeCounter
integer :: neededHalos
Expand All @@ -131,7 +131,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
integer, pointer :: config_n_bcl_iter_mid, config_n_bcl_iter_beg, config_n_bcl_iter_end
integer, pointer :: config_n_ts_iter, config_btr_subcycle_loop_factor
integer, pointer :: config_n_btr_cor_iter, config_num_halos
logical, pointer :: config_use_standardGM
logical, pointer :: config_use_GM,config_use_Redi
integer, pointer :: config_reset_debugTracers_top_nLayers

logical, pointer :: config_use_freq_filtered_thickness, config_btr_solve_SSH2, config_filter_btr_mode
logical, pointer :: config_vel_correction, config_prescribe_velocity, config_prescribe_thickness
Expand All @@ -142,6 +143,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
logical, pointer :: config_use_tracerGroup
logical, pointer :: config_compute_active_tracer_budgets
logical, pointer :: config_use_tidal_potential_forcing
logical, pointer :: config_reset_debugTracers_near_surface

character (len=StrKIND), pointer :: config_land_ice_flux_mode

Expand All @@ -164,6 +166,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

real (kind=RKIND), dimension(:), pointer :: dcEdge, fEdge, bottomDepth, refBottomDepthTopOfCell
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge

! State Array Pointers
Expand Down Expand Up @@ -196,7 +199,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
real (kind=RKIND), dimension(:), pointer :: gradSSHX, gradSSHY, gradSSHZ
real (kind=RKIND), dimension(:), pointer :: gradSSHZonal, gradSSHMeridional
real (kind=RKIND), dimension(:,:), pointer :: surfaceVelocity, SSHGradient
real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta
real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta

! Diagnostics Field Pointers
type (field2DReal), pointer :: normalizedRelativeVorticityEdgeField, divergenceField, relativeVorticityField
Expand Down Expand Up @@ -225,6 +228,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
integer :: threadNum

integer :: temp_mask
real (kind=RKIND) :: tracer2_value, lat

call mpas_timer_start("se timestep")

Expand Down Expand Up @@ -256,13 +260,16 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
call mpas_pool_get_config(domain % configs, 'config_prescribe_velocity', config_prescribe_velocity)
call mpas_pool_get_config(domain % configs, 'config_prescribe_thickness', config_prescribe_thickness)

call mpas_pool_get_config(domain % configs, 'config_use_standardGM', config_use_standardGM)
call mpas_pool_get_config(domain % configs, 'config_use_GM', config_use_GM)
call mpas_pool_get_config(domain % configs, 'config_use_Redi', config_use_Redi)
call mpas_pool_get_config(domain % configs, 'config_use_cvmix_kpp', config_use_cvmix_kpp)
call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode)

call mpas_pool_get_config(domain % configs, 'config_num_halos', config_num_halos)

call mpas_pool_get_config(domain % configs, 'config_compute_active_tracer_budgets', config_compute_active_tracer_budgets)
call mpas_pool_get_config(domain % configs, 'config_reset_debugTracers_near_surface', config_reset_debugTracers_near_surface)
call mpas_pool_get_config(domain % configs, 'config_reset_debugTracers_top_nLayers', config_reset_debugTracers_top_nLayers)
allocate(n_bcl_iter(config_n_ts_iter))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Expand Down Expand Up @@ -1028,7 +1035,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

! Subtract tidal potential from ssh, if needed
! Subtract the tidal potential from the current and new subcycle ssh and store and a work arrays.
! Then point sshSubcycleCur and ssh SubcycleNew to the work arrays so the tidal potential terms
! Then point sshSubcycleCur and ssh SubcycleNew to the work arrays so the tidal potential terms
! are included in the grad operator inside the edge loop.
if (config_use_tidal_potential_forcing) then
call mpas_pool_get_array(forcingPool,'sshSubcycleCurWithTides', sshSubcycleCurWithTides)
Expand Down Expand Up @@ -1397,8 +1404,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! + normalVelocityCorrection
! This is u 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_standardGM) on adding normalGMBolusVelocity
! I think it is not needed because normalGMBolusVelocity=0 when GM not on.
!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) &
*( normalBarotropicVelocityNew(iEdge) + normalBaroclinicVelocityNew(k,iEdge) &
Expand Down Expand Up @@ -1740,6 +1747,58 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp end do
end if
! Reset debugTracers to fixed value at the surface
if ( trim(groupItr % memberName) == 'debugTracers' ) then
call mpas_pool_get_array(meshPool, 'latCell', latCell)
call mpas_pool_get_array(meshPool, 'lonCell', lonCell)
if (config_reset_debugTracers_near_surface) then
!$omp do schedule(runtime)
do iCell = 1, nCells
! Reset tracer1 to 2 in top n layers
do k = 1, config_reset_debugTracers_top_nLayers
tracersGroupNew(1,k,iCell) = 2.0_RKIND
end do
! Reset tracer2 to 2 in top n layers
! in zonal bands, and 1 outside
lat = latCell(iCell)*180./3.1415
if ( lat>-60.0.and.lat<-55.0 &
.or.lat>-40.0.and.lat<-35.0 &
.or.lat>- 2.5.and.lat< 2.5 &
.or.lat> 35.0.and.lat< 40.0 &
.or.lat> 55.0.and.lat< 60.0 ) then
do k = 1, config_reset_debugTracers_top_nLayers
tracersGroupNew(2,k,iCell) = 2.0_RKIND
end do
else
do k = 1, config_reset_debugTracers_top_nLayers
tracersGroupNew(2,k,iCell) = 1.0_RKIND
end do
end if
! Reset tracer3 to 2 in top n layers
! in zonal bands, and 1 outside
lat = latCell(iCell)*180./3.1415
if ( lat>-55.0.and.lat<-50.0 &
.or.lat>-35.0.and.lat<-30.0 &
.or.lat>-15.0.and.lat<-10.0 &
.or.lat> 10.0.and.lat< 15.0 &
.or.lat> 30.0.and.lat< 35.0 &
.or.lat> 50.0.and.lat< 55.0 ) then
do k = 1, config_reset_debugTracers_top_nLayers
tracersGroupNew(3,k,iCell) = 2.0_RKIND
end do
else
do k = 1, config_reset_debugTracers_top_nLayers
tracersGroupNew(3,k,iCell) = 1.0_RKIND
end do
end if
end do
!$omp end do
end if
end if
end if
end if
end do
Expand Down Expand Up @@ -1823,8 +1882,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool)
! Compute normalGMBolusVelocity; it will be added to the baroclinic modes in Stage 2 above.
if (config_use_standardGM) then
call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
if (config_use_GM.or.config_use_Redi) then
call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &
meshPool, scratchPool, timeLevelIn=2)
end if
call ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPool, scratchPool, err, 2)
Expand Down Expand Up @@ -1921,8 +1981,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
call ocn_effective_density_in_land_ice_update(meshPool, forcingPool, statePool, scratchPool, err)

! Compute normalGMBolusVelocity; it will be added to normalVelocity in Stage 2 of the next cycle.
if (config_use_standardGM) then
call ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
if (config_use_GM.or.config_use_Redi) then
call ocn_gm_compute_Bolus_velocity(statePool, diagnosticsPool, &
meshPool, scratchPool, timeLevelIn=2)
end if

call mpas_timer_start('se final mpas reconstruct', .false.)
Expand Down Expand Up @@ -1951,7 +2012,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

call ocn_time_average_coupled_accumulate(diagnosticsPool, statePool, forcingPool, 2)

if (config_use_standardGM) then
if (config_use_GM) then
call ocn_reconstruct_gm_vectors(diagnosticsPool, meshPool)
end if

Expand Down
Loading

0 comments on commit 3e4d773

Please sign in to comment.