From 4cb296c4003014fe57d6d00f86868a78a532fc95 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Tue, 25 Jul 2023 16:11:33 +0000 Subject: [PATCH] Modification of edge mask computation when l_fixed_area=T in horizontal remapping (#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 --- .../cicedyn/dynamics/ice_transport_driver.F90 | 12 +- .../cicedyn/dynamics/ice_transport_remap.F90 | 113 ++++++++---------- cicecore/cicedyn/general/ice_init.F90 | 34 ++++-- doc/source/science_guide/sg_horiztrans.rst | 2 +- 4 files changed, 78 insertions(+), 83 deletions(-) diff --git a/cicecore/cicedyn/dynamics/ice_transport_driver.F90 b/cicecore/cicedyn/dynamics/ice_transport_driver.F90 index 30fe546e0..4f9d84d98 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_driver.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_driver.F90 @@ -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! @@ -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 @@ -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 @@ -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 !------------------------------------------------------------------- diff --git a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 index 286a51711..eb0dd17cf 100644 --- a/cicecore/cicedyn/dynamics/ice_transport_remap.F90 +++ b/cicecore/cicedyn/dynamics/ice_transport_remap.F90 @@ -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 @@ -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. !======================================================================= @@ -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 !======================================================================= @@ -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, & @@ -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) @@ -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 @@ -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), & @@ -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 @@ -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 @@ -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. @@ -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 diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index 2c8b1db3b..3b8d83d1f 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -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, & @@ -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), & @@ -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, & @@ -2770,7 +2771,6 @@ subroutine init_state endif - !----------------------------------------------------------------- ! compute aggregate ice state and open water area !----------------------------------------------------------------- @@ -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, & @@ -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) @@ -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 diff --git a/doc/source/science_guide/sg_horiztrans.rst b/doc/source/science_guide/sg_horiztrans.rst index 7862b5689..10b668755 100644 --- a/doc/source/science_guide/sg_horiztrans.rst +++ b/doc/source/science_guide/sg_horiztrans.rst @@ -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: