Skip to content

Commit

Permalink
Merge pull request mom-ocean#1495 from marshallward/ale_sponge_halo_fix
Browse files Browse the repository at this point in the history 
 ALE sponge mask_z halo fixes
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA authored Sep 13, 2021
2 parents 1623085 + e05ceb2 commit 87aad26
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 53 deletions.
46 changes: 18 additions & 28 deletions src/framework/MOM_horizontal_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,12 +263,16 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
real, intent(in) :: conversion !< Conversion factor for tracer.
integer, intent(in) :: recnum !< Record number of tracer to be read.
type(ocean_grid_type), intent(inout) :: G !< Grid object
real, allocatable, dimension(:,:,:) :: tr_z !< pointer to allocatable tracer array on local
real, allocatable, dimension(:,:,:), intent(out) :: tr_z
!< pointer to allocatable tracer array on local
!! model grid and input-file vertical levels.
real, allocatable, dimension(:,:,:) :: mask_z !< pointer to allocatable tracer mask array on
real, allocatable, dimension(:,:,:), intent(out) :: mask_z
!< pointer to allocatable tracer mask array on
!! local model grid and input-file vertical levels.
real, allocatable, dimension(:) :: z_in !< Cell grid values for input data.
real, allocatable, dimension(:) :: z_edges_in !< Cell grid edge values for input data.
real, allocatable, dimension(:), intent(out) :: z_in
!< Cell grid values for input data.
real, allocatable, dimension(:), intent(out) :: z_edges_in
!< Cell grid edge values for input data.
real, intent(out) :: missing_value !< The missing value in the returned array.
logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
Expand Down Expand Up @@ -329,10 +333,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
is_ongrid = .false.
if (present(ongrid)) is_ongrid = ongrid

if (allocated(tr_z)) deallocate(tr_z)
if (allocated(mask_z)) deallocate(mask_z)
if (allocated(z_edges_in)) deallocate(z_edges_in)

PI_180 = atan(1.0)/45.

! Open NetCDF file and if present, extract data and spatial coordinate information
Expand Down Expand Up @@ -383,13 +383,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion,
rcode = NF90_GET_ATT(ncid, varid, "scale_factor", scale_factor)
if (rcode /= 0) scale_factor = 1.0

if (allocated(lon_in)) deallocate(lon_in)
if (allocated(lat_in)) deallocate(lat_in)
if (allocated(z_in)) deallocate(z_in)
if (allocated(z_edges_in)) deallocate(z_edges_in)
if (allocated(tr_z)) deallocate(tr_z)
if (allocated(mask_z)) deallocate(mask_z)

allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1))
allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))

Expand Down Expand Up @@ -620,12 +613,16 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
type(time_type), intent(in) :: Time !< A FMS time type
real, intent(in) :: conversion !< Conversion factor for tracer.
type(ocean_grid_type), intent(inout) :: G !< Grid object
real, allocatable, dimension(:,:,:) :: tr_z !< pointer to allocatable tracer array on local
real, allocatable, dimension(:,:,:), intent(out) :: tr_z
!< pointer to allocatable tracer array on local
!! model grid and native vertical levels.
real, allocatable, dimension(:,:,:) :: mask_z !< pointer to allocatable tracer mask array on
real, allocatable, dimension(:,:,:), intent(out) :: mask_z
!< pointer to allocatable tracer mask array on
!! local model grid and native vertical levels.
real, allocatable, dimension(:) :: z_in !< Cell grid values for input data.
real, allocatable, dimension(:) :: z_edges_in !< Cell grid edge values for input data. (Intent out)
real, allocatable, dimension(:), intent(out) :: z_in
!< Cell grid values for input data.
real, allocatable, dimension(:), intent(out) :: z_edges_in
!< Cell grid edge values for input data.
real, intent(out) :: missing_value !< The missing value in the returned array.
logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction
logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid
Expand All @@ -651,8 +648,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
integer :: i,j,k
integer, dimension(4) :: start, count, dims, dim_id
real, dimension(:,:), allocatable :: x_in, y_in
real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file
real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole
real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file
real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole
real :: max_lat, min_lat, pole, max_depth, npole
real :: roundoff ! The magnitude of roundoff, usually ~2e-16.
logical :: add_np
Expand Down Expand Up @@ -698,12 +695,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
call cpu_clock_begin(id_clock_read)

call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value)
if (allocated(lon_in)) deallocate(lon_in)
if (allocated(lat_in)) deallocate(lat_in)
if (allocated(z_in)) deallocate(z_in)
if (allocated(z_edges_in)) deallocate(z_edges_in)
if (allocated(tr_z)) deallocate(tr_z)
if (allocated(mask_z)) deallocate(mask_z)

id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3)

Expand Down Expand Up @@ -900,7 +891,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t
enddo
enddo
endif

end subroutine horiz_interp_and_extrap_tracer_fms_id

!> Create a 2d-mesh of grid coordinates from 1-d arrays.
Expand Down
47 changes: 22 additions & 25 deletions src/parameterizations/vertical/MOM_ALE_sponge.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -863,8 +863,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid
real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2]
real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields
real, allocatable, dimension(:,:,:) :: sp_val_u ! A temporary array for fields
real, allocatable, dimension(:,:,:) :: sp_val_v ! A temporary array for fields
real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts
real, allocatable, dimension(:,:,:) :: mask_u ! A temporary array for field mask at u pts
real, allocatable, dimension(:,:,:) :: mask_v ! A temporary array for field mask at v pts
Expand All @@ -883,6 +881,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m].
real :: sp_val_u ! Interpolation of sp_val to u-points
real :: sp_val_v ! Interpolation of sp_val to v-points
integer :: nPoints

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Expand All @@ -903,8 +903,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
call MOM_error(FATAL,"apply_ALE_sponge: No time information provided")
do m=1,CS%fldno
nz_data = CS%Ref_val(m)%nz_data
allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)); sp_val(:,:,:) = 0.0
allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)); mask_z(:,:,:) = 0.0
call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, &
z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, &
Expand Down Expand Up @@ -991,24 +989,20 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
call MOM_error(FATAL,"apply_ALE_sponge: No time information provided")

nz_data = CS%Ref_val_u%nz_data
allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
allocate(sp_val_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
allocate(mask_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
sp_val(:,:,:) = 0.0
sp_val_u(:,:,:) = 0.0
mask_u(:,:,:) = 0.0
mask_z(:,:,:) = 0.0
! Interpolate from the external horizontal grid and in time
call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, 1.0, G, sp_val, mask_z, z_in, &
z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,&
answers_2018=CS%hor_regrid_answers_2018)

! Initialize mask_z halos to zero before pass_var, in case of no update
mask_z(G%isc-1, G%jsc:G%jec, :) = 0.
mask_z(G%iec+1, G%jsc:G%jec, :) = 0.
call pass_var(sp_val, G%Domain)
call pass_var(mask_z, G%Domain)

allocate(mask_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
do j=G%jsc,G%jec; do I=G%iscB,G%iecB
sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data))
mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data))
enddo ; enddo

Expand All @@ -1018,7 +1012,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
! Therefore we use c as per C code and increment the index where necessary.
i = CS%col_i_u(c) ; j = CS%col_j_u(c)
if (mask_u(i,j,1) == 1.0) then
CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data)
do k=1,nz_data
sp_val_u = 0.5 * (sp_val(i,j,k) + sp_val(i+1,j,k))
CS%Ref_val_u%p(k,c) = sp_val_u
enddo
else
CS%Ref_val_u%p(1:nz_data,c) = 0.0
endif
Expand All @@ -1039,24 +1036,21 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) )
CS%Ref_val_u%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data)
enddo
deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc)
deallocate(sp_val, mask_u, mask_z, hsrc)
nz_data = CS%Ref_val_v%nz_data
allocate(sp_val( G%isd:G%ied,G%jsd:G%jed,1:nz_data))
allocate(sp_val_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
allocate(mask_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
sp_val(:,:,:) = 0.0
sp_val_v(:,:,:) = 0.0
mask_z(:,:,:) = 0.0
! Interpolate from the external horizontal grid and in time
call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, 1.0, G, sp_val, mask_z, z_in, &
z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,&
answers_2018=CS%hor_regrid_answers_2018)
! Initialize mask_z halos to zero before pass_var, in case of no update
mask_z(G%isc:G%iec, G%jsc-1, :) = 0.
mask_z(G%isc:G%iec, G%jec+1, :) = 0.
call pass_var(sp_val, G%Domain)
call pass_var(mask_z, G%Domain)

allocate(mask_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
do J=G%jscB,G%jecB; do i=G%isc,G%iec
sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data))
mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data))
enddo ; enddo
!call pass_var(mask_z,G%Domain)
Expand All @@ -1066,7 +1060,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
! Therefore we use c as per C code and increment the index where necessary.
i = CS%col_i_v(c) ; j = CS%col_j_v(c)
if (mask_v(i,j,1) == 1.0) then
CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data)
do k=1,nz_data
sp_val_v = 0.5 * (sp_val(i,j,k) + sp_val(i,j+1,k))
CS%Ref_val_v%p(k,c) = sp_val_v
enddo
else
CS%Ref_val_v%p(1:nz_data,c) = 0.0
endif
Expand All @@ -1087,7 +1084,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) )
CS%Ref_val_v%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data)
enddo
deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc)
deallocate(sp_val, mask_v, mask_z, hsrc)
endif

call pass_var(h,G%Domain)
Expand Down

0 comments on commit 87aad26

Please sign in to comment.