Skip to content

Commit

Permalink
Fix Redi resolution taper
Browse files Browse the repository at this point in the history
  • Loading branch information
mark-petersen committed Apr 21, 2022
1 parent 88fd9c0 commit 0996cfa
Show file tree
Hide file tree
Showing 2 changed files with 11 additions and 8 deletions.
9 changes: 5 additions & 4 deletions components/mpas-ocean/src/shared/mpas_ocn_gm.F
Original file line number Diff line number Diff line change
Expand Up @@ -940,7 +940,7 @@ subroutine ocn_GM_init(domain, err)!{{{

integer, intent(out) :: err !< Output: error flag

real(kind=RKIND) :: sphereRadius
real(kind=RKIND) :: sphereRadius, coef

type(block_type), pointer :: block
type(mpas_pool_type), pointer :: meshPool
Expand Down Expand Up @@ -1075,7 +1075,7 @@ subroutine ocn_GM_init(domain, err)!{{{
call mpas_dmpar_finalize(domain%dminfo)
end if

! Add horizontal taper
! Initialize horizontal taper
if (config_GM_horizontal_taper == 'none' .or. &
config_GM_horizontal_taper == 'RossbyRadius') then
! For 'RossbyRadius', the taper is recomputed at every time step.
Expand All @@ -1088,6 +1088,8 @@ subroutine ocn_GM_init(domain, err)!{{{
!$omp end do
!$omp end parallel
else if (config_GM_horizontal_taper == 'ramp') then
coef = 1.0_RKIND &
/(config_GM_horizontal_ramp_max - config_GM_horizontal_ramp_min)
!$omp parallel
!$omp do schedule(runtime)
do iEdge = 1, nEdges
Expand All @@ -1096,8 +1098,7 @@ subroutine ocn_GM_init(domain, err)!{{{
else if (dcEdge(iEdge) >= config_GM_horizontal_ramp_max) then
gmHorizontalTaper(iEdge) = 1.0_RKIND
else
gmHorizontalTaper(iEdge) = 1.0_RKIND &
/(config_GM_horizontal_ramp_max - config_GM_horizontal_ramp_min) &
gmHorizontalTaper(iEdge) = coef &
*(dcEdge(iEdge) - config_GM_horizontal_ramp_min)
end if
end do
Expand Down
10 changes: 6 additions & 4 deletions components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
Original file line number Diff line number Diff line change
Expand Up @@ -506,6 +506,7 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{
real(kind=RKIND), dimension(:), pointer :: RediKappa, RediKappaData
real(kind=RKIND), dimension(:), pointer :: RediHorizontalTaper
real(kind=RKIND), dimension(:), pointer :: dcEdge
real(kind=RKIND) :: coef
integer, dimension(:,:), pointer :: cellsOnEdge

integer :: k, iEdge
Expand Down Expand Up @@ -533,7 +534,7 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)

! Resolution taper
! Initialize horizonal taper
if (config_Redi_horizontal_taper == 'none' .or. &
config_Redi_horizontal_taper == 'RossbyRadius') then
! For 'RossbyRadius', the taper is recomputed at every time step in
Expand All @@ -546,6 +547,8 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{
!$omp end do
!$omp end parallel
else if (config_Redi_horizontal_taper == 'ramp') then
coef = 1.0_RKIND &
/(config_Redi_horizontal_ramp_max - config_Redi_horizontal_ramp_min)
!$omp parallel
!$omp do schedule(runtime)
do iEdge=1,nEdges
Expand All @@ -554,9 +557,8 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{
else if (dcEdge(iEdge) >= config_Redi_horizontal_ramp_max) then
RediHorizontalTaper(iEdge) = 1.0_RKIND
else
RediHorizontalTaper(iEdge) = RediHorizontalTaper(iEdge) &
*(dcEdge(iEdge) - config_Redi_horizontal_ramp_min) &
/(config_Redi_horizontal_ramp_max - config_Redi_horizontal_ramp_min)
RediHorizontalTaper(iEdge) = coef &
*(dcEdge(iEdge) - config_Redi_horizontal_ramp_min)
end if
end do
!$omp end do
Expand Down

0 comments on commit 0996cfa

Please sign in to comment.