Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ALE sponge mask_z halo fixes #1495

Merged
merged 4 commits into from
Sep 13, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 18 additions & 28 deletions src/framework/MOM_horizontal_regridding.F90
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
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