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.

Sign up for GitHub

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

Jump to bottom

Modification of edge mask computation when l_fixed_area=T in horizontal remapping #833

Merged
merged 9 commits into from
Jul 25, 2023
Merged

Modification of edge mask computation when l_fixed_area=T in horizontal remapping #833

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
86d3409
Use same method whether l_fixed_area=T or F to compute masks for edge…
JFLemieux73 Jun 2, 2023
76c9bb4
Corrected typo in comment
JFLemieux73 Jun 5, 2023
21ad9c0
Cosmetic (indentation) change in ice_transport_remap.F90
JFLemieux73 Jun 7, 2023
68f768c
Set l_fixed_area value depending of grid type
JFLemieux73 Jun 12, 2023
3d26d4c
Modifs to the doc for l_fixed_area
JFLemieux73 Jun 12, 2023
f93d322
Use umask for uvel,vvel initialization for boxslotcyl and change grid…
JFLemieux73 Jun 19, 2023
387c7df
Temporary changes before next PR: l_fixed_area=F for B and C grid
JFLemieux73 Jul 24, 2023
816dda0
Temporary changes before next PR: remove paragraph in the doc
JFLemieux73 Jul 24, 2023
84149cc
Small modifs: l_fixed_area and grid_ice are defined in module ice_tra…
JFLemieux73 Jul 24, 2023
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
12 changes: 4 additions & 8 deletions cicecore/cicedyn/dynamics/ice_transport_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,6 @@ module ice_transport_driver
! 'upwind' => 1st order donor cell scheme
! 'remap' => remapping scheme

logical, parameter :: &
l_fixed_area = .false. ! if true, prescribe area flux across each edge

! NOTE: For remapping, hice and hsno are considered tracers.
! ntrace is not equal to ntrcr!

Expand Down Expand Up @@ -81,6 +78,7 @@ subroutine init_transport
use ice_state, only: trcr_depend
use ice_timers, only: ice_timer_start, ice_timer_stop, timer_advect
use ice_transport_remap, only: init_remap
use ice_grid, only: grid_ice

integer (kind=int_kind) :: &
k, nt, nt1 ! tracer indices
Expand Down Expand Up @@ -236,7 +234,7 @@ subroutine init_transport
endif ! master_task
1000 format (1x,a,2x,i6,2x,i6,2x,i4,4x,l4)

if (trim(advection)=='remap') call init_remap ! grid quantities
if (trim(advection)=='remap') call init_remap ! grid quantities

call ice_timer_stop(timer_advect) ! advection

Expand Down Expand Up @@ -545,19 +543,17 @@ subroutine transport_remap (dt)
call horizontal_remap (dt, ntrace, &
uvel (:,:,:), vvel (:,:,:), &
aim (:,:,:,:), trm(:,:,:,:,:), &
l_fixed_area, &
tracer_type, depend, &
has_dependents, integral_order, &
l_dp_midpt, grid_ice, &
l_dp_midpt, &
uvelE (:,:,:), vvelN (:,:,:))
else
call horizontal_remap (dt, ntrace, &
uvel (:,:,:), vvel (:,:,:), &
aim (:,:,:,:), trm(:,:,:,:,:), &
l_fixed_area, &
tracer_type, depend, &
has_dependents, integral_order, &
l_dp_midpt, grid_ice)
l_dp_midpt)
endif

!-------------------------------------------------------------------
Expand Down
113 changes: 52 additions & 61 deletions cicecore/cicedyn/dynamics/ice_transport_remap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module ice_transport_remap
use ice_exit, only: abort_ice
use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
use icepack_intfc, only: icepack_query_parameters
use ice_grid, only : grid_ice

implicit none
private
Expand All @@ -57,6 +58,11 @@ module ice_transport_remap
p5625m = -9._dbl_kind/16._dbl_kind ,&
p52083 = 25._dbl_kind/48._dbl_kind

logical :: &
l_fixed_area ! if true, prescribe area flux across each edge
! if false, area flux is determined internally
! and is passed out

logical (kind=log_kind), parameter :: bugcheck = .false.

!=======================================================================
Expand Down Expand Up @@ -293,6 +299,29 @@ subroutine init_remap
enddo
!$OMP END PARALLEL DO

!-------------------------------------------------------------------
! Set logical l_fixed_area depending of the grid type.
!
! If l_fixed_area is true, the area of each departure region is
! computed in advance (e.g., by taking the divergence of the
! velocity field and passed to locate_triangles. The departure
! regions are adjusted to obtain the desired area.
! If false, edgearea is computed in locate_triangles and passed out.
!
! l_fixed_area = .false. has been the default approach in CICE. It is
! used like this for the B-grid. However, idealized tests with the
! C-grid have shown that l_fixed_area = .false. leads to a checkerboard
! pattern in prognostic fields (e.g. aice). Using l_fixed_area = .true.
! eliminates the checkerboard pattern in C-grid simulations.
!
!-------------------------------------------------------------------

if (grid_ice == 'CD' .or. grid_ice == 'C') then
l_fixed_area = .false. !jlem temporary
else
l_fixed_area = .false.
endif

end subroutine init_remap

!=======================================================================
Expand All @@ -316,11 +345,10 @@ end subroutine init_remap
subroutine horizontal_remap (dt, ntrace, &
uvel, vvel, &
mm, tm, &
l_fixed_area, &
tracer_type, depend, &
has_dependents, &
integral_order, &
l_dp_midpt, grid_ice, &
l_dp_midpt, &
uvelE, vvelN)

use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
Expand Down Expand Up @@ -353,21 +381,6 @@ subroutine horizontal_remap (dt, ntrace, &
real (kind=dbl_kind), intent(inout), dimension (nx_block,ny_block,ntrace,ncat,max_blocks) :: &
tm ! mean tracer values in each grid cell

character (len=char_len_long), intent(in) :: &
grid_ice ! ice grid, B, C, etc

!-------------------------------------------------------------------
! If l_fixed_area is true, the area of each departure region is
! computed in advance (e.g., by taking the divergence of the
! velocity field and passed to locate_triangles. The departure
! regions are adjusted to obtain the desired area.
! If false, edgearea is computed in locate_triangles and passed out.
!-------------------------------------------------------------------

logical, intent(in) :: &
l_fixed_area ! if true, edgearea_e and edgearea_n are prescribed
! if false, edgearea is computed here and passed out

integer (kind=int_kind), dimension (ntrace), intent(in) :: &
tracer_type , & ! = 1, 2, or 3 (see comments above)
depend ! tracer dependencies (see above)
Expand Down Expand Up @@ -716,8 +729,7 @@ subroutine horizontal_remap (dt, ntrace, &
dxu (:,:,iblk), dyu(:,:,iblk), &
xp (:,:,:,:), yp (:,:,:,:), &
iflux, jflux, &
triarea, &
l_fixed_area, edgearea_e(:,:))
triarea, edgearea_e(:,:))

!-------------------------------------------------------------------
! Given triangle vertices, compute coordinates of triangle points
Expand Down Expand Up @@ -776,8 +788,7 @@ subroutine horizontal_remap (dt, ntrace, &
dxu (:,:,iblk), dyu (:,:,iblk), &
xp (:,:,:,:), yp(:,:,:,:), &
iflux, jflux, &
triarea, &
l_fixed_area, edgearea_n(:,:))
triarea, edgearea_n(:,:))

call triangle_coordinates (nx_block, ny_block, &
integral_order, icellsng(:,iblk), &
Expand Down Expand Up @@ -1696,8 +1707,7 @@ subroutine locate_triangles (nx_block, ny_block, &
dxu, dyu, &
xp, yp, &
iflux, jflux, &
triarea, &
l_fixed_area, edgearea)
triarea, edgearea)

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
Expand Down Expand Up @@ -1730,12 +1740,6 @@ subroutine locate_triangles (nx_block, ny_block, &
indxi , & ! compressed index in i-direction
indxj ! compressed index in j-direction

logical, intent(in) :: &
l_fixed_area ! if true, the area of each departure region is
! passed in as edgearea
! if false, edgearea if determined internally
! and is passed out

real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
edgearea ! area of departure region for each edge
! edgearea > 0 for eastward/northward flow
Expand Down Expand Up @@ -1838,7 +1842,7 @@ subroutine locate_triangles (nx_block, ny_block, &
! BL | BC | BR (bottom left, center, right)
! | |
!
! and the transport is across the edge between cells TC and TB.
! and the transport is across the edge between cells TC and BC.
!
! Departure points are scaled to a local coordinate system
! whose origin is at the midpoint of the edge.
Expand Down Expand Up @@ -1951,45 +1955,32 @@ subroutine locate_triangles (nx_block, ny_block, &
! Compute mask for edges with nonzero departure areas
!-------------------------------------------------------------------

if (l_fixed_area) then
icellsd = 0
icellsd = 0
if (trim(edge) == 'north') then
do j = jb, je
do i = ib, ie
if (edgearea(i,j) /= c0) then
if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
else
icellsd = 0
if (trim(edge) == 'north') then
do j = jb, je
do i = ib, ie
if (dpx(i-1,j)/=c0 .or. dpy(i-1,j)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
else ! east edge
do j = jb, je
do i = ib, ie
if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
endif ! edge = north/east
endif ! l_fixed_area
else ! east edge
do j = jb, je
do i = ib, ie
if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 &
.or. &
dpx(i,j)/=c0 .or. dpy(i,j)/=c0) then
icellsd = icellsd + 1
indxid(icellsd) = i
indxjd(icellsd) = j
endif
enddo
enddo
endif ! edge = north/east

!-------------------------------------------------------------------
! Scale the departure points
Expand Down
34 changes: 21 additions & 13 deletions cicecore/cicedyn/general/ice_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2532,7 +2532,7 @@ subroutine init_state
use ice_domain, only: nblocks, blocks_ice, halo_info
use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero, nfsd
use ice_flux, only: sst, Tf, Tair, salinz, Tmltz
use ice_grid, only: tmask, ULON, TLAT, grid_ice, grid_average_X2Y
use ice_grid, only: tmask, umask, ULON, TLAT, grid_ice, grid_average_X2Y
use ice_boundary, only: ice_HaloUpdate
use ice_constants, only: field_loc_Nface, field_loc_Eface, field_type_scalar
use ice_state, only: trcr_depend, aicen, trcrn, vicen, vsnon, &
Expand Down Expand Up @@ -2720,6 +2720,7 @@ subroutine init_state
ilo, ihi, jlo, jhi, &
iglob, jglob, &
ice_ic, tmask(:,:, iblk), &
umask(:,:, iblk), &
ULON (:,:, iblk), &
TLAT (:,:, iblk), &
Tair (:,:, iblk), sst (:,:, iblk), &
Expand All @@ -2742,10 +2743,10 @@ subroutine init_state

if (grid_ice == 'CD' .or. grid_ice == 'C') then

call grid_average_X2Y('S',uvel,'U',uvelN,'N')
call grid_average_X2Y('S',vvel,'U',vvelN,'N')
call grid_average_X2Y('S',uvel,'U',uvelE,'E')
call grid_average_X2Y('S',vvel,'U',vvelE,'E')
call grid_average_X2Y('A',uvel,'U',uvelN,'N')
call grid_average_X2Y('A',vvel,'U',vvelN,'N')
call grid_average_X2Y('A',uvel,'U',uvelE,'E')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this change in the averaging type needed?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is just for the initialization of boxslotcyl test. uvelE and vvelN are defined from the uvel, vvel values. uvel(78,79) is initialized to 0 because of the land mask (no slip) while uvel(78,78) is initialized to 2.27. As uvelE is interpolated from uvel, uvelE(78,79) should be 2.27/2 =1.135. This is obtained by changing the averaging type.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, thank you.

call grid_average_X2Y('A',vvel,'U',vvelE,'E')

! Halo update on North, East faces
call ice_HaloUpdate(uvelN, halo_info, &
Expand All @@ -2760,7 +2761,6 @@ subroutine init_state

endif


!-----------------------------------------------------------------
! compute aggregate ice state and open water area
!-----------------------------------------------------------------
Expand Down Expand Up @@ -2819,8 +2819,9 @@ subroutine set_state_var (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
iglob, jglob, &
ice_ic, tmask, &
ULON, &
TLAT, &
umask, &
ULON, &
TLAT, &
Tair, sst, &
Tf, &
salinz, Tmltz, &
Expand All @@ -2845,7 +2846,8 @@ subroutine set_state_var (nx_block, ny_block, &
ice_ic ! method of ice cover initialization

logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: &
tmask ! true for ice/ocean cells
tmask , & ! true for ice/ocean cells
umask ! for U points

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
ULON , & ! longitude of velocity pts (radians)
Expand Down Expand Up @@ -3293,13 +3295,19 @@ subroutine set_state_var (nx_block, ny_block, &
domain_length = dxrect*cm_to_m*nx_global
period = c12*secday ! 12 days rotational period
max_vel = pi*domain_length/period

do j = 1, ny_block
do i = 1, nx_block

uvel(i,j) = c2*max_vel*(real(jglob(j), kind=dbl_kind) - p5) &
/ real(ny_global - 1, kind=dbl_kind) - max_vel
vvel(i,j) = -c2*max_vel*(real(iglob(i), kind=dbl_kind) - p5) &
/ real(nx_global - 1, kind=dbl_kind) + max_vel
if (umask(i,j)) then
uvel(i,j) = c2*max_vel*(real(jglob(j), kind=dbl_kind) - p5) &
/ real(ny_global - 1, kind=dbl_kind) - max_vel
vvel(i,j) = -c2*max_vel*(real(iglob(i), kind=dbl_kind) - p5) &
/ real(nx_global - 1, kind=dbl_kind) + max_vel
else
uvel(i,j) = c0
vvel(i,j) = c0
endif
enddo ! j
enddo ! i
else
Expand Down
2 changes: 1 addition & 1 deletion doc/source/science_guide/sg_horiztrans.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ remapping scheme of :cite:`Dukowicz00` as modified for sea ice by

- The upwind scheme uses velocity points at the East and North face (i.e. :math:`uvelE=u` at the E point and :math:`vvelN=v` at the N point) of a T gridcell. As such, the prognostic C grid velocity components (:math:`uvelE` and :math:`vvelN`) can be passed directly to the upwind transport scheme. If the upwind scheme is used with the B grid, the B grid velocities, :math:`uvelU` and :math:`vvelU` (respectively :math:`u` and :math:`v` at the U point) are interpolated to the E and N points first. (Note however that the upwind scheme does not transport all potentially available tracers.)

- The remapping scheme uses :math:`uvelU` and :math:`vvelU` if l_fixed_area is false and :math:`uvelE` and :math:`vvelN` if l_fixed_area is true. l_fixed_area is hardcoded to false by default and further described below. As such, the B grid velocities (:math:`uvelU` and :math:`vvelU`) are used directly in the remapping scheme, while the C grid velocities (:math:`uvelE` and :math:`vvelN`) are interpolated to U points first. If l_fixed_area is changed to true, then the reverse is true. The C grid velocities are used directly and the B grid velocities are interpolated.
- Remapping is naturally a B-grid transport scheme as the corner (U point) velocity components :math:`uvelU` and :math:`vvelU` are used to calculate departure points. Nevertheless, the remapping scheme can also be used with the C grid by first interpolating :math:`uvelE` and :math:`vvelN` to the U points.

The remapping scheme has several desirable features:

Expand Down
Loading