Skip to content

Commit

Permalink
Modification of edge mask computation when l_fixed_area=T in horizont…
Browse files Browse the repository at this point in the history 
…al remapping (CICE-Consortium#833)

* Use same method whether l_fixed_area=T or F to compute masks for edge fluxes

* Corrected typo in comment

* Cosmetic (indentation) change in ice_transport_remap.F90

* Set l_fixed_area value depending of grid type

* Modifs to the doc for l_fixed_area

* Use umask for uvel,vvel initialization for boxslotcyl and change grid avg type from S to A in init_state

* Temporary changes before next PR: l_fixed_area=F for B and C grid

* Temporary changes before next PR: remove paragraph in the doc

* Small modifs: l_fixed_area and grid_ice are defined in module ice_transport_remap
  • Loading branch information
@JFLemieux73
JFLemieux73 authored Jul 25, 2023
1 parent 9f42a62 commit 4cb296c
Show file tree
Hide file tree
Showing 4 changed files with 78 additions and 83 deletions.
12 changes: 4 additions & 8 deletions cicecore/cicedyn/dynamics/ice_transport_driver.F90
View file
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
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
Original file line number Diff line number Diff line change
Expand Up @@ -2542,7 +2542,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 @@ -2730,6 +2730,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 @@ -2752,10 +2753,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')
call grid_average_X2Y('A',vvel,'U',vvelE,'E')

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

endif


!-----------------------------------------------------------------
! compute aggregate ice state and open water area
!-----------------------------------------------------------------
Expand Down Expand Up @@ -2829,8 +2829,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 @@ -2855,7 +2856,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 @@ -3303,13 +3305,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
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

0 comments on commit 4cb296c

Please sign in to comment.