From 67d16df77245228aecf42523cf2d04b0a1092c0f Mon Sep 17 00:00:00 2001 From: dengwirda Date: Tue, 25 Apr 2023 00:32:09 -0600 Subject: [PATCH 01/19] Update SAL and AMHarmonicAnalysis to respect initial SSH offsets - Include pressure forcing term in SAL computation to account for under isc and atm pressure cases. - Update SAL coastal smoothing factor to mask shallow cells where wet-dry effects may occur out of SAL computation. - Compute harmonic analysis AM wrt. SSH-SSH_init to account for initially non-zero SSH distributions. --- components/mpas-ocean/src/Registry.xml | 14 ++-- .../Registry_harmonic_analysis.xml | 3 + .../mpas_ocn_harmonic_analysis.F | 36 ++++++--- .../mpas_ocn_time_integration_si.F | 18 ++--- .../mpas_ocn_time_integration_split.F | 20 ++--- .../src/mode_init/mpas_ocn_init_test_sht.F | 17 ++-- .../shared/mpas_ocn_diagnostics_variables.F | 12 +-- .../mpas-ocean/src/shared/mpas_ocn_tendency.F | 5 +- .../mpas_ocn_vel_self_attraction_loading.F | 78 ++++++++++--------- .../src/shared/mpas_ocn_vel_tidal_potential.F | 22 +++--- 10 files changed, 127 insertions(+), 98 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index bb399577eee3..90a218151a2f 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -829,8 +829,8 @@ description="Controls if self-attraction and loading is applied to ssh" possible_values=".true. or .false." /> - - + - - diff --git a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml index 9337cc1e6e23..6fcfc5b0c560 100644 --- a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml +++ b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml @@ -80,6 +80,9 @@ + diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F index 6f53e885f784..ee1495a3e59a 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F @@ -117,11 +117,12 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{ type (dm_info) :: dminfo type (block_type), pointer :: block - type (mpas_pool_type), pointer :: harmonicAnalysisAMPool + type (mpas_pool_type), pointer :: statePool, harmonicAnalysisAMPool - integer, pointer :: nAnalysisConstituents + integer, pointer :: nAnalysisConstituents, nCellsSolve real (kind=RKIND), pointer :: harmonicAnalysisStart real (kind=RKIND), pointer :: harmonicAnalysisEnd + real (kind=RKIND), dimension(:), pointer :: ssh, sshInit real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentAmplitude @@ -135,7 +136,7 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{ integer, dimension(:), allocatable :: tidalConstituentType type (MPAS_Time_Type) :: refTime - integer :: iCon + integer :: iCon, iCell err = 0 @@ -145,8 +146,13 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{ block => domain % blocklist do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'state', statePool) call mpas_pool_get_subpool(block % structs, 'harmonicAnalysisAM', harmonicAnalysisAMPool) + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + + call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit) call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents) call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix) call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector) @@ -214,11 +220,14 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{ tidalConstituentType, & err) - if (.not. config_do_restart) then - leastSquaresRHSVector = 0.0_RKIND - leastSquaresLHSMatrix = 0.0_RKIND - end if + if (.not. config_do_restart) then + leastSquaresRHSVector = 0.0_RKIND + leastSquaresLHSMatrix = 0.0_RKIND + end if + do iCell = 1,nCellSSolve + sshInit(iCell) = ssh(iCell) + end do do iCon = 1,nAnalysisConstituents call mpas_log_write('Constituent '//constituentList(iCon)%constituent) @@ -291,7 +300,8 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{ integer, pointer :: nCellsSolve integer, pointer :: nAnalysisConstituents - real (kind=RKIND), dimension(:), pointer :: ssh + real (kind=RKIND), dimension(:), pointer :: ssh, sshInit + real (kind=RKIND), dimension(:), allocatable :: sshDiff real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector real (kind=RKIND), dimension(:), pointer :: analysisConstituentFrequency @@ -334,6 +344,7 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{ call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit) call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents) call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix) call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector) @@ -355,6 +366,12 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{ ! update if within harmonic analysis period if (daysSinceStartOfSim .le. harmonicAnalysisEnd) then + allocate(sshDiff(nCellsSolve)) + + do iCell = 1,nCellsSolve + sshDiff(iCell) = ssh(iCell) - sshInit(iCell) + end do + CALL update_least_squares_LHS_matrix(nAnalysisConstituents, & time, & analysisConstituentFrequency, & @@ -364,8 +381,9 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{ time, & nCellsSolve, & analysisConstituentFrequency, & - ssh, & + sshDiff, & leastSquaresRHSVector) + deallocate(sshDiff) call mpas_log_write('harmonicAnalysisAM update') end if 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 743c250ed306..ae43f3a084bf 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 @@ -155,7 +155,7 @@ module ocn_time_integration_si SIvec_zh1_field,SIvec_wh0_field,SIvec_rh1_field ! field pointer for halo update - integer :: ssh_sal_on + integer :: pgf_sal_on #ifdef USE_GPU_AWARE_MPI ! SI - halo exch related variables ---------------------------------- @@ -715,7 +715,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ endif if (config_use_tidal_potential_forcing) then - !$acc enter data copyin(ssh_sal) + !$acc enter data copyin(pgf_sal) !$acc enter data create(sshSubcycleCurWithTides,tidalPotentialEta) endif #endif @@ -1160,7 +1160,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$acc parallel loop gang & !$acc present(sshSubcycleCurWithTides,sshSubcycleCur, & - !$acc tidalPotentialEta,ssh_sal,sshSubcycleCur) + !$acc tidalPotentialEta,pgf_sal,sshSubcycleCur) #else !$omp parallel !$omp do schedule(runtime) @@ -1169,8 +1169,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ sshSubcycleCurWithTides(iCell) = & sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & - ssh_sal_on * ssh_sal(iCell) - & - (1 - ssh_sal_on) * self_attraction_and_loading_beta* & + pgf_sal_on * pgf_sal(iCell) - & + (1 - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) end do #ifndef MPAS_OPENACC @@ -2801,7 +2801,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ #ifdef MPAS_OPENACC if (config_use_tidal_potential_forcing) then - !$acc exit data delete(ssh_sal,tidalPotentialEta) + !$acc exit data delete(pgf_sal,tidalPotentialEta) !$acc exit data copyout(sshSubcycleCurWithTides) endif #endif @@ -3501,10 +3501,10 @@ subroutine ocn_time_integration_si_init(domain, dt)!{{{ !-------------------------------------------------------------------- if (config_use_self_attraction_loading) then - ! set ssh_sal_on to 1 - ssh_sal_on = 1 + ! set pgf_sal_on to 1 + pgf_sal_on = 1 else - ssh_sal_on = 0 + pgf_sal_on = 0 endif self_attraction_and_loading_beta = config_self_attraction_and_loading_beta 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 4f24b544033c..d294afdc5bfd 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 @@ -91,7 +91,7 @@ module ocn_time_integration_split numTSIterations, &! number of outer timestep iterations nBtrSubcycles ! number of barotropic subcycles - integer :: ssh_sal_on + integer :: pgf_sal_on integer :: & thickEdgeFluxChoice ! choice of thickness flux type integer, parameter :: & @@ -1020,8 +1020,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshSubcycleCurWithTides(iCell) = & sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & - ssh_sal_on * ssh_sal(iCell) - & - (1 - ssh_sal_on) * self_attraction_and_loading_beta* & + pgf_sal_on * pgf_sal(iCell) - & + (1 - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) end do !$omp end do @@ -1268,15 +1268,15 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshSubcycleCurWithTides(iCell) = & sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & - ssh_sal_on * ssh_sal(iCell) - & - (1 - ssh_sal_on) * self_attraction_and_loading_beta* & + pgf_sal_on * pgf_sal(iCell) - & + (1 - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) sshSubcycleNewWithTides(iCell) = & sshSubcycleNew(iCell) - & tidalPotentialEta(iCell) - & - ssh_sal_on * ssh_sal(iCell) - & - (1 - ssh_sal_on) * self_attraction_and_loading_beta* & + pgf_sal_on * pgf_sal(iCell) - & + (1 - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleNew(iCell) end do !$omp end do @@ -3071,10 +3071,10 @@ subroutine ocn_time_integration_split_init(domain)!{{{ !----------------------------------------------------------------- if (config_use_self_attraction_loading) then - ! set ssh_sal_on to 1 - ssh_sal_on = 1 + ! set pgf_sal_on to 1 + pgf_sal_on = 1 else - ssh_sal_on = 0 + pgf_sal_on = 0 endif self_attraction_and_loading_beta = config_self_attraction_and_loading_beta diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F index 33eaf9af0127..ae3450695435 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F @@ -89,7 +89,7 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{ integer, pointer :: nCellsSolve real (kind=RKIND), dimension(:), pointer :: lonCell, latCell real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell - real (kind=RKIND), dimension(:), pointer :: ssh, ssh_sal, ssh_sal_grad + real (kind=RKIND), dimension(:), pointer :: ssh, surfacePressure, pgf_sal, pgf_sal_grad real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge real (kind=RKIND), pointer :: config_test_sht_cosine_bell_lat_center, & config_test_sht_cosine_bell_lon_center, & @@ -188,9 +188,11 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{ call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(statePool, 'ssh', ssh) - call mpas_pool_get_array(diagnosticsPool, 'ssh_sal', ssh_sal) - call mpas_pool_get_array(diagnosticsPool, 'ssh_sal_grad', ssh_sal_grad) + call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', surfacePressure) + call mpas_pool_get_array(diagnosticsPool, 'pgf_sal', pgf_sal) + call mpas_pool_get_array(diagnosticsPool, 'pgf_sal_grad', pgf_sal_grad) + surfacePressure(:) = 0.0_RKIND ! Setup function to approximate if (config_test_sht_function_option == 1) then @@ -283,22 +285,23 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{ call ocn_vel_self_attraction_loading_init(domain,err) do i = 1,config_test_sht_n_iterations call mpas_log_write('Iteration = $i',intArgs=(/i/)) - call ocn_compute_self_attraction_loading(domain, forcingPool, domain % dminfo, ssh, err) + call ocn_compute_self_attraction_loading(domain, forcingPool, domain % dminfo, & + ssh, surfacePressure, err) enddo ! Compute gradient do iEdge = 1,nEdgesArray(1) cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - ssh_sal_grad(iEdge) = (ssh_sal(cell2) - ssh_sal(cell1))/dcEdge(iEdge) + pgf_sal_grad(iEdge) = (pgf_sal(cell2) - pgf_sal(cell1))/dcEdge(iEdge) enddo ! Compute RMS Error error_local_area = 0.0_RKIND error_local = 0.0_RKIND do iCell = 1,nCellsSolve - error_local_area = error_local_area + areaCell(iCell)*(ssh(iCell)-ssh_sal(iCell))**2 - error_local = error_local + (ssh(iCell)-ssh_sal(iCell))**2 + error_local_area = error_local_area + areaCell(iCell)*(ssh(iCell)-pgf_sal(iCell))**2 + error_local = error_local + (ssh(iCell)-pgf_sal(iCell))**2 enddo call mpas_dmpar_sum_real(domain % dminfo,error_local_area,error_global_area) call mpas_dmpar_sum_real(domain % dminfo,error_local,error_global) 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 567c3efc183c..4b8c317c397f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -72,7 +72,7 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:), pointer :: gradSSH real (kind=RKIND), dimension(:), pointer :: pressureAdjustedSSH - real (kind=RKIND), dimension(:), pointer :: ssh_sal + real (kind=RKIND), dimension(:), pointer :: pgf_sal real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceValue real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceLayerValue @@ -375,8 +375,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e if (config_use_tidal_potential_forcing) then call mpas_pool_get_array(diagnosticsPool, & - 'ssh_sal', & - ssh_sal) + 'pgf_sal', & + pgf_sal) end if if ( trim(config_ocean_run_mode) == 'init' .or. & @@ -691,7 +691,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc enter data create(topographic_wave_drag) end if if (config_use_self_attraction_loading) then - !$acc enter data copyin(ssh_sal) + !$acc enter data copyin(pgf_sal) end if if ( trim(config_ocean_run_mode) == 'init' ) then !$acc enter data create( & @@ -938,7 +938,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc exit data delete(topographic_wave_drag) end if if (config_use_self_attraction_loading) then - !$acc exit data delete(ssh_sal) + !$acc exit data delete(pgf_sal) end if if ( trim(config_ocean_run_mode) == 'init' ) then !$acc exit data delete( & @@ -1275,7 +1275,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ nullify(salinitySurfaceRestoringTendency) end if if ( config_use_self_attraction_loading ) then - nullify(ssh_sal) + nullify(pgf_sal) end if if ( trim(config_time_integrator) == 'split_implicit' ) then nullify(SIvec_r0, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index 91f182db9993..a16e52f14499 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -426,12 +426,13 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & thermExpCoeff, salineContractCoeff, & tendVel, err) - ! Compute ssh additions due to self-attraction and loading + ! Compute pgf additions due to self-attraction and loading if(mpas_is_alarm_ringing(domain % clock, 'salComputeAlarm', ierr=err)) then #ifdef MPAS_DEBUG call mpas_log_write(' Computing SAL') #endif - call ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, err) + call ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, & + ssh, surfacePressure, err) call mpas_reset_clock_alarm(domain % clock, 'salComputeAlarm', ierr=err) endif diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F index 6b91daac4d37..4f9be1aef45c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F @@ -128,7 +128,8 @@ module ocn_vel_self_attraction_loading ! !----------------------------------------------------------------------- - subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, err)!{{{ + subroutine ocn_compute_self_attraction_loading(domain, forcingPool, & + dminfo, ssh, surfacePressure, err)!{{{ !----------------------------------------------------------------- ! @@ -137,6 +138,7 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, !----------------------------------------------------------------- type (dm_info), intent(in) :: dminfo real (kind=RKIND), dimension(:), intent(in) :: ssh + real (kind=RKIND), dimension(:), intent(in) :: surfacePressure !----------------------------------------------------------------- ! @@ -161,7 +163,7 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, !----------------------------------------------------------------- integer :: iCell - + call mpas_timer_start('SAL Calculation') err = 0 @@ -171,7 +173,8 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, endif do iCell = 1,nCellsAll - sshSmoothed(iCell) = coastalSmoothingFactor(iCell)*ssh(iCell) + sshSmoothed(iCell) = coastalSmoothingFactor(iCell) * ( & + ssh(iCell) + rho0gInv * surfacePressure(iCell) ) enddo if (.not. config_use_parallel_self_attraction_loading) then @@ -184,9 +187,9 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, endif - ! ssh_sal currently computed on host but must be transferred + ! pgf_sal currently computed on host but must be transferred ! to device for later tendency calculations - !$acc update device(ssh_sal) + !$acc update device(pgf_sal) call mpas_timer_stop('SAL Calculation') @@ -236,7 +239,7 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ !----------------------------------------------------------------- ! SAL variables - real (kind=RKIND), dimension(:), pointer :: ssh_gg + real (kind=RKIND), dimension(:), pointer :: pgf_gg ! MPI Variables integer :: iCell, ilm, curProc @@ -253,7 +256,7 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ allocate(globalArray(nCellsGlobal), gatheredArray(nCellsGlobal)) endif - ! Gather only the nCellsOwned from ssh (does not include Halos) + ! Gather only the nCellsOwned from pgf (does not include Halos) call MPI_GATHERV(sshSmoothed, nCellsOwned, MPI_DOUBLE, gatheredArray, nCellsPerProc, & nCellsDisplacement, MPI_DOUBLE, 0, dminfo % comm, err_tmp) err = ior(err, err_tmp) @@ -262,20 +265,20 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ ! Perform SAL calculation only on process 0 if (curProc.eq.0) then call mpas_timer_start('Serial SAL: Forward') - allocate(ssh_gg(nGrid)) - ssh_gg(:) = 0.0_dp + allocate(pgf_gg(nGrid)) + pgf_gg(:) = 0.0_dp Slm(:)=(0.0_dp, 0.0_dp) - ! Rearrange ssh into CellID order + ! Rearrange into CellID order do iCell = 1,nCellsGlobal globalArray(indexToCellIDGathered(iCell)) = gatheredArray(iCell) enddo - ! Interpolate ssh onto Gaussian Grid - call interpolate( toColValues, toRowValues, toSValues, globalArray, ssh_gg) + ! Interpolate pgf onto Gaussian Grid + call interpolate( toColValues, toRowValues, toSValues, globalArray, pgf_gg) ! Perform spherical harmonic transform #ifdef USE_SHTNS - call Spat_to_SH(shtns_c, ssh_gg, Slm) + call Spat_to_SH(shtns_c, pgf_gg, Slm) call mpas_timer_stop('Serial SAL: Forward') call mpas_timer_start('Serial SAL: Inverse') ! Multiply each harmonic coefficient by the scaling factor @@ -283,10 +286,10 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ Slm(ilm) = Slm(ilm) * LoveScaling(ilm) enddo ! Perform inverse spherical harmonic transform - call SH_to_spat(shtns_c, Slm, ssh_gg) + call SH_to_spat(shtns_c, Slm, pgf_gg) #endif ! Interpolate back to MPAS mesh - call interpolate( fromColValues, fromRowValues, fromSValues, ssh_gg, globalArray) + call interpolate( fromColValues, fromRowValues, fromSValues, pgf_gg, globalArray) ! Rearrange back to index order do iCell = 1,nCellsGlobal @@ -300,18 +303,18 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ call mpas_timer_stop('Serial SAL: Inverse') endif - ! Scatter back to ssh_sal + ! Scatter back to pgf_sal call mpas_timer_start('Serial SAL: Scatter') call MPI_SCATTERV(gatheredArray, nCellsPerProc, nCellsDisplacement, MPI_DOUBLE, & - ssh_sal, nCellsAll, MPI_DOUBLE, 0, dminfo % comm, err_tmp) + pgf_sal, nCellsAll, MPI_DOUBLE, 0, dminfo % comm, err_tmp) err = ior(err, err_tmp) if (curProc.eq.0) then - deallocate(globalArray, gatheredArray, ssh_gg) + deallocate(globalArray, gatheredArray, pgf_gg) endif ! Perform Halo exchange update - call mpas_dmpar_field_halo_exch(domain,'ssh_sal') + call mpas_dmpar_field_halo_exch(domain,'pgf_sal') call mpas_timer_stop('Serial SAL: Scatter') end subroutine self_attraction_loading_serial!}}} @@ -514,7 +517,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ !!!!!!!!!!!!!!!!!!!! do iCell = 1,nCellsAll - ssh_sal(iCell) = 0.0_RKIND + pgf_sal(iCell) = 0.0_RKIND enddo do m = 0,nOrder @@ -539,7 +542,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Sum together product of spherical harmonic functions and coefficients do iCell = endIdx, startIdx, -1 - ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmnm2(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) + pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmnm2(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) enddo !------------ @@ -553,7 +556,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Sum together product of spherical harmonic functions and coefficients do iCell = endIdx, startIdx, -1 - ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmnm1(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) + pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmnm1(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) enddo endif @@ -576,7 +579,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Sum together product of spherical harmonic functions and coefficients do iCell = endIdx, startIdx, -1 - ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmn(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) + pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmn(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1)) enddo enddo ! n loop @@ -609,6 +612,7 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{ ! Config variables type (block_type), pointer :: block_ptr type (mpas_pool_type), pointer :: forcingPool + type (mpas_pool_type), pointer :: statePool ! NetCDF and weights file variables integer :: toNcId, toNsDimId, toRowId, toColId, toSId @@ -641,8 +645,8 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{ integer :: nCellBlock integer :: err_tmp - real(kind=RKIND) :: d - real(kind=RKIND) :: trans_start, trans_width + real(kind=RKIND), dimension(:), pointer :: ssh + real(kind=RKIND) :: depth, depth_tol call mpas_timer_start('SAL Init') @@ -665,24 +669,24 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{ block_ptr => domain % blocklist call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) + call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) call mpas_pool_get_array(forcingPool, 'coastalSmoothingFactor', & coastalSmoothingFactor) + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) - ! Set up coastal ssh smoothing + ! Set up coastal smoothing allocate(sshSmoothed(nCellsAll)) - if (config_self_attraction_loading_smoothing_width > 2.0_RKIND) then - trans_width = config_self_attraction_loading_smoothing_width*1000.0_RKIND - trans_start = 0.0_RKIND - do iCell = 1,nCellsAll - d = distanceToCoast(iCell) - coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh((d - trans_start - 0.5_RKIND * trans_width) / (0.2_RKIND *trans_width)) + 1.0_RKIND) - enddo - else - coastalSmoothingFactor(:) = 1.0_RKIND - endif - + ! Set-up smoothing mask for shallow regions where SAL should not + ! be computed across e.g. wet-dry cells + depth_tol = config_self_attraction_loading_depth_cutoff + do iCell = 1,nCellsAll + depth = ssh(iCell) + bottomDepth(iCell) + coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh( & + (depth - 0.5_RKIND*depth_tol) / (0.2_RKIND*depth_tol)) + 1.0_RKIND) + enddo + ! Setup Alarm for SAL alarmTime = mpas_get_clock_time(domain % clock, MPAS_START_TIME, ierr=err) call mpas_set_timeInterval(alarmTimeStep, timeString = & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F index 092c52c54d09..211538ac0a73 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F @@ -95,8 +95,8 @@ module ocn_vel_tidal_potential type(char_array), dimension(37) :: constituentList public :: char_array - ! Flag for ssh_sal - real(kind=RKIND) :: ssh_sal_on + ! Flag for pgf_sal + real(kind=RKIND) :: pgf_sal_on !*********************************************************************** @@ -162,7 +162,7 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{ #ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, & - !$acc tidalPotEta, ssh, surfacePressure, tend, ssh_sal) & + !$acc tidalPotEta, ssh, surfacePressure, tend, pgf_sal) & !$acc private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) #else !$omp parallel do schedule(runtime) & @@ -176,8 +176,8 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{ potentialGrad = - gravity*invdcEdge* & ( tidalPotEta(cell2) - tidalPotEta(cell1) & - + (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * ( rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) & - + ssh_sal_on * (ssh_sal(cell2) - ssh_sal(cell1)) ) + + (1.0_RKIND - pgf_sal_on) * betaSelfAttrLoad * ( rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) & + + pgf_sal_on * (pgf_sal(cell2) - pgf_sal(cell1)) ) do k=1,kMax tend(k,iEdge) = tend(k,iEdge) - & @@ -193,7 +193,7 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{ #ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, & - !$acc tidalPotEta, ssh, surfacePressure, tend, ssh_sal) & + !$acc tidalPotEta, ssh, surfacePressure, tend, pgf_sal) & !$acc private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) #else !$omp parallel do schedule(runtime) & @@ -207,9 +207,9 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{ potentialGrad = - gravity*invdcEdge* & ( tidalPotEta(cell2) - tidalPotEta(cell1) & - + (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad *( (ssh(cell2) - ssh(cell1)) & - + rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) & - + ssh_sal_on * (ssh_sal(cell2) - ssh_sal(cell1)) ) + + (1.0_RKIND - pgf_sal_on) * betaSelfAttrLoad *( (ssh(cell2) - ssh(cell1)) & + + rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) & + + pgf_sal_on * (pgf_sal(cell2) - pgf_sal(cell1)) ) do k=1,kMax tend(k,iEdge) = tend(k,iEdge) - & @@ -366,9 +366,9 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ tidalPotentialOff = .false. if (config_use_self_attraction_loading) then - ssh_sal_on = 1.0_RKIND + pgf_sal_on = 1.0_RKIND else - ssh_sal_on = 0.0_RKIND + pgf_sal_on = 0.0_RKIND endif betaSelfAttrLoad = config_self_attraction_and_loading_beta From 14e19d35207ee2a330a9b984e692deabce933de5 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Tue, 25 Apr 2023 00:58:45 -0600 Subject: [PATCH 02/19] Update SAL and AMHarmonicAnalysis to respect initial SSH offsets - Add missing rho0gInv definition. --- .../src/shared/mpas_ocn_vel_self_attraction_loading.F | 3 +++ 1 file changed, 3 insertions(+) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F index 4f9be1aef45c..e425b49110d9 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F @@ -163,6 +163,7 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, & !----------------------------------------------------------------- integer :: iCell + real (kind=RKIND) :: rho0gInv call mpas_timer_start('SAL Calculation') @@ -172,6 +173,8 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, & return endif + rho0gInv = 1.0_RKIND / rho_sw / gravity + do iCell = 1,nCellsAll sshSmoothed(iCell) = coastalSmoothingFactor(iCell) * ( & ssh(iCell) + rho0gInv * surfacePressure(iCell) ) From 992db9bf3f87f29529c98c2d3cd4b6824d38ab30 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Tue, 25 Apr 2023 13:50:27 -0600 Subject: [PATCH 03/19] Update SAL and AMHarmonicAnalysis to respect initial SSH offsets - Clean up descriptions, comments and use real-valued switches. --- .../analysis_members/Registry_harmonic_analysis.xml | 3 ++- .../src/mode_forward/mpas_ocn_time_integration_si.F | 8 ++++---- .../mode_forward/mpas_ocn_time_integration_split.F | 12 ++++++------ .../shared/mpas_ocn_vel_self_attraction_loading.F | 2 ++ 4 files changed, 14 insertions(+), 11 deletions(-) diff --git a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml index 6fcfc5b0c560..a43f3c0536d3 100644 --- a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml +++ b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml @@ -81,7 +81,7 @@ description="Frequency of each constituent" /> + 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 ae43f3a084bf..9529fbdc6ea1 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 @@ -155,7 +155,7 @@ module ocn_time_integration_si SIvec_zh1_field,SIvec_wh0_field,SIvec_rh1_field ! field pointer for halo update - integer :: pgf_sal_on + real(kind=RKIND) :: pgf_sal_on #ifdef USE_GPU_AWARE_MPI ! SI - halo exch related variables ---------------------------------- @@ -1170,7 +1170,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & pgf_sal_on * pgf_sal(iCell) - & - (1 - pgf_sal_on) * self_attraction_and_loading_beta* & + (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) end do #ifndef MPAS_OPENACC @@ -3502,9 +3502,9 @@ subroutine ocn_time_integration_si_init(domain, dt)!{{{ if (config_use_self_attraction_loading) then ! set pgf_sal_on to 1 - pgf_sal_on = 1 + pgf_sal_on = 1.0_RKIND else - pgf_sal_on = 0 + pgf_sal_on = 0.0_RKIND endif self_attraction_and_loading_beta = config_self_attraction_and_loading_beta 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 d294afdc5bfd..653d88198131 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 @@ -91,12 +91,12 @@ module ocn_time_integration_split numTSIterations, &! number of outer timestep iterations nBtrSubcycles ! number of barotropic subcycles - integer :: pgf_sal_on integer :: & thickEdgeFluxChoice ! choice of thickness flux type integer, parameter :: & thickEdgeFluxCenter = 1, &! use mean thickness of cell neighbors thickEdgeFluxUpwind = 2 ! use upwind cell thickness at edge + real (kind=RKIND) :: pgf_sal_on real (kind=RKIND) :: & useVelocityCorrection,&! mask for velocity correction @@ -1021,7 +1021,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & pgf_sal_on * pgf_sal(iCell) - & - (1 - pgf_sal_on) * self_attraction_and_loading_beta* & + (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) end do !$omp end do @@ -1269,14 +1269,14 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshSubcycleCur(iCell) - & tidalPotentialEta(iCell) - & pgf_sal_on * pgf_sal(iCell) - & - (1 - pgf_sal_on) * self_attraction_and_loading_beta* & + (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleCur(iCell) sshSubcycleNewWithTides(iCell) = & sshSubcycleNew(iCell) - & tidalPotentialEta(iCell) - & pgf_sal_on * pgf_sal(iCell) - & - (1 - pgf_sal_on) * self_attraction_and_loading_beta* & + (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* & sshSubcycleNew(iCell) end do !$omp end do @@ -3072,9 +3072,9 @@ subroutine ocn_time_integration_split_init(domain)!{{{ !----------------------------------------------------------------- if (config_use_self_attraction_loading) then ! set pgf_sal_on to 1 - pgf_sal_on = 1 + pgf_sal_on = 1.0_RKIND else - pgf_sal_on = 0 + pgf_sal_on = 0.0_RKIND endif self_attraction_and_loading_beta = config_self_attraction_and_loading_beta diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F index e425b49110d9..c34115fb3ff6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F @@ -176,6 +176,8 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, & rho0gInv = 1.0_RKIND / rho_sw / gravity do iCell = 1,nCellsAll + ! an adjusted ssh that accounts for surface pressure, + ! with SAL computed wrt the full btr bottom pressure sshSmoothed(iCell) = coastalSmoothingFactor(iCell) * ( & ssh(iCell) + rho0gInv * surfacePressure(iCell) ) enddo From 3c32bd47a1b2f936451b7ad72d834fd7701e9a4c Mon Sep 17 00:00:00 2001 From: dengwirda Date: Tue, 25 Apr 2023 14:03:26 -0600 Subject: [PATCH 04/19] Update SAL and AMHarmonicAnalysis to respect initial SSH offsets - Change sshSmoothed to btrPressure in SAL routines, since we're now working with surface pressures too. --- .../mpas_ocn_vel_self_attraction_loading.F | 35 +++++++++---------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F index c34115fb3ff6..90d1a4e4630d 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F @@ -103,7 +103,7 @@ module ocn_vel_self_attraction_loading real(kind=RKIND), dimension(:), allocatable :: Snm, SnmRe, SnmIm real(kind=RKIND), dimension(:), allocatable :: Snm_local, SnmRe_local, SnmIm_local real(kind=RKIND), dimension(:,:), allocatable :: Snm_local_reproSum, SnmRe_local_reproSum, SnmIm_local_reproSum - real(kind=RKIND), dimension(:), allocatable :: sshSmoothed + real(kind=RKIND), dimension(:), allocatable :: btrPressure real(kind=RKIND), dimension(:), pointer :: coastalSmoothingFactor integer, dimension(:,:), allocatable :: blockIdxForward integer, dimension(:,:), allocatable :: blockIdxInverse @@ -176,9 +176,8 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, & rho0gInv = 1.0_RKIND / rho_sw / gravity do iCell = 1,nCellsAll - ! an adjusted ssh that accounts for surface pressure, - ! with SAL computed wrt the full btr bottom pressure - sshSmoothed(iCell) = coastalSmoothingFactor(iCell) * ( & + ! compute SAL wrt. the full barotropic bottom pressure + btrPressure(iCell) = coastalSmoothingFactor(iCell) * ( & ssh(iCell) + rho0gInv * surfacePressure(iCell) ) enddo @@ -262,7 +261,7 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{ endif ! Gather only the nCellsOwned from pgf (does not include Halos) - call MPI_GATHERV(sshSmoothed, nCellsOwned, MPI_DOUBLE, gatheredArray, nCellsPerProc, & + call MPI_GATHERV(btrPressure, nCellsOwned, MPI_DOUBLE, gatheredArray, nCellsPerProc, & nCellsDisplacement, MPI_DOUBLE, 0, dminfo % comm, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop('Serial SAL: Gather') @@ -402,13 +401,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Compute local integral contribution if (config_parallel_self_attraction_loading_bfb) then do iCell = endIdx, startIdx, -1 - SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1) enddo else do iCell = endIdx, startIdx, -1 - SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1) enddo endif @@ -424,13 +423,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Compute local integral contribution if (config_parallel_self_attraction_loading_bfb) then do iCell = endIdx, startIdx, -1 - SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1) enddo else do iCell = endIdx, startIdx, -1 - SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1) enddo endif @@ -455,13 +454,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{ ! Compute local integral contribution if (config_parallel_self_attraction_loading_bfb) then do iCell = endIdx, startIdx, -1 - SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1) enddo else do iCell = endIdx, startIdx, -1 - SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1) - SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1) + SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1) + SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1) enddo endif @@ -681,7 +680,7 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{ call mpas_pool_get_array(statePool, 'ssh', ssh, 1) ! Set up coastal smoothing - allocate(sshSmoothed(nCellsAll)) + allocate(btrPressure(nCellsAll)) ! Set-up smoothing mask for shallow regions where SAL should not ! be computed across e.g. wet-dry cells From 75ebe3df42925306e7b99d555f927abe8f6a174d Mon Sep 17 00:00:00 2001 From: Kristin Barton Date: Fri, 18 Nov 2022 12:19:47 -0800 Subject: [PATCH 05/19] Added two new topographic wave drag schemes. --- components/mpas-ocean/src/Registry.xml | 30 +++- .../shared/mpas_ocn_diagnostics_variables.F | 32 ++++- .../src/shared/mpas_ocn_vel_forcing.F | 2 +- ...as_ocn_vel_forcing_topographic_wave_drag.F | 128 +++++++++++++++++- 4 files changed, 178 insertions(+), 14 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 90a218151a2f..c1f5c23669ff 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1108,7 +1108,11 @@ description="If true, topographic wave drag is used on the momentum equation." possible_values=".true. or .false." /> - + @@ -2262,7 +2266,11 @@ filename_template="topographic_wave_drag.nc" runtime_format="single_file" mode="forward;analysis"> - + + + + + @@ -3233,10 +3241,26 @@ - + + + + \brief Computes tendency term from topographic wave drag !> \author Nairita Pal @@ -75,7 +77,7 @@ module ocn_vel_forcing_topographic_wave_drag ! !----------------------------------------------------------------------- - subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{ + subroutine ocn_vel_forcing_topographic_wave_drag_scalar_tend(normalVelocity, & !{{{ tend, err) !----------------------------------------------------------------- @@ -117,6 +119,11 @@ subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{ call mpas_timer_start('vel topographic wave drag') + ! JSL: C = pi std(H)^2 N_B / L + ! ZAE: C = gamma H (grad(H))^2 N_B bar(N) / (8 * pi^2 * omega) + ! N(z) = N_0 exp(z/1300) + ! Drag: rho_0 C u + !$omp do schedule(runtime) private(k) do iEdge = 1, nEdgesOwned k = maxLevelEdgeTop(iEdge) @@ -126,7 +133,7 @@ subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{ ! topographic_wave_drag = 1/rinv if (k>0) then tend(k,iEdge) = tend(k,iEdge) - topographicWaveDragCoeff & - * topographic_wave_drag(iEdge) * normalVelocity(k,iEdge) + * topographicWaveDrag(iEdge) * normalVelocity(k,iEdge) endif enddo !$omp end do @@ -135,7 +142,85 @@ subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{ !-------------------------------------------------------------------- - end subroutine ocn_vel_forcing_topographic_wave_drag_tend!}}} + end subroutine ocn_vel_forcing_topographic_wave_drag_scalar_tend!}}} + +!*********************************************************************** +! +! routine ocn_vel_forcing_topographic_wave_drag_tensor_tend +! +!> \brief Computes tendency term from topographic wave drag +!> \author Kristin Barton +!> \date 1 Nov 2022 +!> \details +!> This routine computes the topographic wave drag tendency for momentum +!> based on current state using the tensor formulation +! +!*********************************************************************** + + subroutine ocn_vel_forcing_topographic_wave_drag_tensor_tend(normalVelocity, & !{{{ + tend, err) + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalVelocity !< Input: velocity + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k + + err = 0 + if ( .not. topographicWaveDragOn ) return + + call mpas_timer_start('vel topographic wave drag') + + ! JSL: C = pi std(H)^2 N_B / L + ! ZAE: C = gamma H (grad(H))^2 N_B bar(N) / (8 * pi^2 * omega) + ! N(z) = N_0 exp(z/1300) + ! Drag: rho_0 C u + + !$omp do schedule(runtime) private(k) + do iEdge = 1, nEdgesOwned + k = maxLevelEdgeTop(iEdge) + ! topographic wave drag term: + ! du/dt = ... 1/rinv * u + ! rinv is the e-folding time use in HyCOM. Here + ! topographic_wave_drag = 1/rinv + if (k>0) then + endif + enddo + !$omp end do + + call mpas_timer_stop('vel topographic wave drag') + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_forcing_topographic_wave_drag_tensor_tend!}}} !*********************************************************************** ! @@ -169,6 +254,39 @@ subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{ if (config_use_topographic_wave_drag) then topographicWaveDragOn = .true. topographicWaveDragCoeff = config_topographic_wave_drag_coeff + if (config_topographic_wave_drag_scheme.EQ."JSL") then + call mpas_log_write("") + call mpas_log_write(" Topographic wave drag scheme is: Jayne and St. Laurent") + call mpas_log_write("") + ! C = pi std(H)^2 Nb / L, L=10,000 + topographicWaveDrag = pii*(bathy_stddev**2)*topo_buoyancy_N1B/10000.0_RKIND + else if (config_topographic_wave_drag_scheme.EQ."ZAE") then + call mpas_log_write("") + call mpas_log_write(" Topographic wave drag scheme is: Zaron and Egbert") + call mpas_log_write("") + ! C = G*H*(grad H)^2 * Nb * Nv / (8 pi^2 omega), take M2 omega = 1.405189e-4_RKIND +! Nv = (5.24e-3_RKIND)*1300*(1 - exp(-bottomDepth/1300))/bottomDepth +! Nb = (5.24e-3_RKIND)*exp(-bottomDepth/1300) + topographicWaveDrag = bottomDepth*(bed_slope_edges**2)*((5.24e-3_RKIND)*1300.0_RKIND*(1 - exp(-bottomDepth/1300.0_RKIND))/bottomDepth)*((5.24e-3_RKIND)*exp(-bottomDepth/1300.0_RKIND))/ (8 * 1.405189e-4_RKIND * pii**2) +! topographicWaveDrag = bottomDepth*(bed_slope_edges**2)*topo_buoyancy_N1B*topo_buoyancy_N1V / (8 * 1.405189e-4_RKIND * pii**2) + else if (config_topographic_wave_drag_scheme.EQ."NYC") then + call mpas_log_write("") + call mpas_log_write(" Topographic wave drag scheme is: Nycander") + call mpas_log_write("") + topographicWaveDrag = topo_rinv + ! Preperation for Nycander + else + call mpas_log_write("") + call mpas_log_write("Invalid parameter for config_topopgraphic_wave_drag_scheme. It must be one of: 'JSL', 'ZAE', or 'NYC'" , MPAS_LOG_CRIT) + ! C = Cit sqrt((Nb^2 - omega^2)(Nm^2 - omega^2))/(4 pi omega) * [ [grad(hl)^2, grad(hl)*grad(hp) ], [ grad(hl) * grad(hp), grad(hp)^2 ] ] + ! Cit = tuning coefficint + ! Nb = bottom frequency + ! Nm = depth-averaged frequency + ! grad(hl) = zonal topo gradient + ! grad(hp) = meridional topo gradient + + call mpas_log_write("") + end if endif if (config_disable_vel_topographic_wave_drag) topographicWaveDragOn = .false. From d02d7405b8776b79813a47984c391c3d307a4496 Mon Sep 17 00:00:00 2001 From: Kristin Barton Date: Mon, 13 Feb 2023 13:30:35 -0800 Subject: [PATCH 06/19] Added local generation formula topograhpic wave drag scheme --- components/mpas-ocean/src/Registry.xml | 47 ++++- .../shared/mpas_ocn_diagnostics_variables.F | 53 +++++- .../src/shared/mpas_ocn_vel_forcing.F | 2 +- ...as_ocn_vel_forcing_topographic_wave_drag.F | 167 ++++++------------ 4 files changed, 153 insertions(+), 116 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index c1f5c23669ff..e92b831fd116 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1109,8 +1109,8 @@ possible_values=".true. or .false." /> + + + + + + + + + @@ -3261,6 +3276,18 @@ description="The integral gradient average per cell" packages="topographicWaveDragPKG" /> + + + + + + + .false. -1.0 +10.0 'mpas_to_grid.nc' 'grid_to_mpas.nc' '0000-00-00_00:30:00' @@ -399,8 +399,11 @@ 1.0e-1 1.0e-3 .false. +'ZAE' 5.0e-4 'centered' +500 +10 0.0 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 196d4abf730c..ae872caf67d9 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1229,9 +1229,9 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml - -Defines region over which ssh is smoothed to zero at coasts for SAL calculation. +Defines depths over which ssh is smoothed to zero for SAL calculation. Valid values: Any positive real number. Default: Defined in namelist_defaults.xml @@ -1780,11 +1780,19 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml + +Which of the following three wave drag schemes to use: JSL: Jayne and St. Laurent; ZAE: Zaron and Egbert; or LGF: Local generation formula + +Valid values: JSL or ZAE or LGF +Default: Defined in namelist_defaults.xml + + Dimensionless topographic wave drag coefficient, $c_{topo}$. -Valid values: any positive real, typically 5.0e-4 +Valid values: any positive real Default: Defined in namelist_defaults.xml @@ -1796,6 +1804,22 @@ Valid values: 'harmonic-mean', 'centered' Default: Defined in namelist_defaults.xml + +Specifies minimum depth at which topographic wave drag is applied + +Valid values: any positive real +Default: Defined in namelist_defaults.xml + + + +Width over which the topographic wave drag is progressively turned off + +Valid values: any positive real +Default: Defined in namelist_defaults.xml + +