From f810375e1b5c297cc14b8e81f65982b822cec170 Mon Sep 17 00:00:00 2001 From: apcraig Date: Wed, 1 Mar 2023 13:48:48 -0700 Subject: [PATCH 01/13] Initial halochk unit test implementation - Add halochk unit test - Add "unknown" and "noupdate" checks to ice_boundary - Remove field_loc_Wface, not used anywhere, not supported - Update cice_decomp.csh script To Do: validate tripole and tripoleT, add unit test to test suite --- .../infrastructure/comm/mpi/ice_boundary.F90 | 191 +++++- .../comm/serial/ice_boundary.F90 | 191 +++++- .../drivers/unittest/halochk/CICE_InitMod.F90 | 472 ++++++++++++++ cicecore/drivers/unittest/halochk/halochk.F90 | 608 ++++++++++++++++++ cicecore/shared/ice_constants.F90 | 3 +- configuration/scripts/Makefile | 6 +- configuration/scripts/cice_decomp.csh | 6 +- configuration/scripts/options/set_env.halochk | 2 + configuration/scripts/options/set_nml.bcyclic | 3 + configuration/scripts/options/set_nml.bopen | 3 + doc/source/user_guide/ug_testing.rst | 1 + 11 files changed, 1478 insertions(+), 8 deletions(-) create mode 100644 cicecore/drivers/unittest/halochk/CICE_InitMod.F90 create mode 100644 cicecore/drivers/unittest/halochk/halochk.F90 create mode 100644 configuration/scripts/options/set_env.halochk create mode 100644 configuration/scripts/options/set_nml.bcyclic create mode 100644 configuration/scripts/options/set_nml.bopen diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 index 68436cd0f..77af8ceee 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 @@ -18,8 +18,10 @@ module ice_boundary use ice_communicate, only: my_task, mpiR4, mpiR8, mpitagHalo use ice_constants, only: field_type_scalar, & field_type_vector, field_type_angle, & + field_type_unknown, field_type_noupdate, & field_loc_center, field_loc_NEcorner, & - field_loc_Nface, field_loc_Eface + field_loc_Nface, field_loc_Eface, & + field_loc_unknown, field_loc_noupdate use ice_global_reductions, only: global_maxval use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -1229,6 +1231,23 @@ subroutine ice_HaloUpdate2DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1653,6 +1672,23 @@ subroutine ice_HaloUpdate2DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2051,6 +2087,23 @@ subroutine ice_HaloUpdate2DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2426,6 +2479,23 @@ subroutine ice_HaloUpdate2DL1(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! copy logical into integer array and call haloupdate on integer array @@ -2519,6 +2589,23 @@ subroutine ice_HaloUpdate3DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2976,6 +3063,23 @@ subroutine ice_HaloUpdate3DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3433,6 +3537,23 @@ subroutine ice_HaloUpdate3DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3890,6 +4011,23 @@ subroutine ice_HaloUpdate4DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -4371,6 +4509,23 @@ subroutine ice_HaloUpdate4DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -4852,6 +5007,23 @@ subroutine ice_HaloUpdate4DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -5319,6 +5491,23 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate_stress)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index 2b81c4441..65da2a6f4 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -17,8 +17,10 @@ module ice_boundary use ice_communicate, only: my_task use ice_constants, only: field_type_scalar, & field_type_vector, field_type_angle, & + field_type_unknown, field_type_noupdate, & field_loc_center, field_loc_NEcorner, & - field_loc_Nface, field_loc_Eface + field_loc_Nface, field_loc_Eface, & + field_loc_unknown, field_loc_noupdate use ice_global_reductions, only: global_maxval use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -699,6 +701,23 @@ subroutine ice_HaloUpdate2DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1011,6 +1030,23 @@ subroutine ice_HaloUpdate2DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1303,6 +1339,23 @@ subroutine ice_HaloUpdate2DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1583,6 +1636,23 @@ subroutine ice_HaloUpdate2DL1(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! copy logical into integer array and call haloupdate on integer array @@ -1662,6 +1732,23 @@ subroutine ice_HaloUpdate3DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1980,6 +2067,23 @@ subroutine ice_HaloUpdate3DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2298,6 +2402,23 @@ subroutine ice_HaloUpdate3DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2616,6 +2737,23 @@ subroutine ice_HaloUpdate4DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2951,6 +3089,23 @@ subroutine ice_HaloUpdate4DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3286,6 +3441,23 @@ subroutine ice_HaloUpdate4DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3610,6 +3782,23 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate_stress)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value diff --git a/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 b/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 new file mode 100644 index 000000000..9ed1c5cbc --- /dev/null +++ b/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 @@ -0,0 +1,472 @@ +!======================================================================= +! +! This module contains the CICE initialization routine that sets model +! parameters and initializes the grid and CICE state variables. +! +! authors Elizabeth C. Hunke, LANL +! William H. Lipscomb, LANL +! Philip W. Jones, LANL +! +! 2006: Converted to free form source (F90) by Elizabeth Hunke +! 2008: E. Hunke moved ESMF code to its own driver + + module CICE_InitMod + + use ice_kinds_mod + use ice_exit, only: abort_ice + use ice_fileunits, only: init_fileunits, nu_diag + use icepack_intfc, only: icepack_aggregate + use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist + use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_configure + use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted + use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, & + icepack_query_tracer_indices, icepack_query_tracer_sizes + + implicit none + private + public :: CICE_Initialize, cice_init + +!======================================================================= + + contains + +!======================================================================= + +! Initialize the basic state, grid and all necessary parameters for +! running the CICE model. Return the initial state in routine +! export state. +! Note: This initialization driver is designed for standalone and +! CESM-coupled applications. For other +! applications (e.g., standalone CAM), this driver would be +! replaced by a different driver that calls subroutine cice_init, +! where most of the work is done. + + subroutine CICE_Initialize + + character(len=*), parameter :: subname='(CICE_Initialize)' + !-------------------------------------------------------------------- + ! model initialization + !-------------------------------------------------------------------- + + call cice_init + + end subroutine CICE_Initialize + +!======================================================================= +! +! Initialize CICE model. + + subroutine cice_init + + use ice_arrays_column, only: hin_max, c_hi_range, alloc_arrays_column + use ice_arrays_column, only: floe_rad_l, floe_rad_c, & + floe_binwidth, c_fsd_range + use ice_state, only: alloc_state + use ice_flux_bgc, only: alloc_flux_bgc + use ice_calendar, only: dt, dt_dyn, istep, istep1, write_ic, & + init_calendar, advance_timestep, calc_timesteps + use ice_communicate, only: init_communicate, my_task, master_task + use ice_diagnostics, only: init_diags + use ice_domain, only: init_domain_blocks + use ice_domain_size, only: ncat, nfsd + use ice_dyn_eap, only: init_eap + use ice_dyn_evp, only: init_evp + use ice_dyn_vp, only: init_vp + use ice_dyn_shared, only: kdyn + use ice_flux, only: init_coupler_flux, init_history_therm, & + init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux + use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, & + get_forcing_atmo, get_forcing_ocn, get_wave_spec + use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & + faero_default, faero_optics, alloc_forcing_bgc, fiso_default + use ice_grid, only: init_grid1, init_grid2, alloc_grid + use ice_history, only: init_hist, accum_hist + use ice_restart_shared, only: restart, runtype + use ice_init, only: input_data, init_state + use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers + use ice_kinds_mod + use ice_restoring, only: ice_HaloRestore_init + use ice_timers, only: timer_total, init_ice_timers, ice_timer_start + use ice_transport_driver, only: init_transport + + logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, & + tr_iso, tr_fsd, wave_spec + character(len=*), parameter :: subname = '(cice_init)' + + call init_communicate ! initial setup for message passing + call init_fileunits ! unit numbers + + ! tcx debug, this will create a different logfile for each pe + ! if (my_task /= master_task) nu_diag = 100+my_task + + call icepack_configure() ! initialize icepack + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + call input_data ! namelist variables + call input_zbgc ! vertical biogeochemistry namelist + call count_tracers ! count tracers + + call init_domain_blocks ! set up block decomposition + call init_grid1 ! domain distribution + call alloc_grid ! allocate grid arrays + call alloc_arrays_column ! allocate column arrays + call alloc_state ! allocate state arrays + call alloc_flux_bgc ! allocate flux_bgc arrays + call alloc_flux ! allocate flux arrays + call init_ice_timers ! initialize all timers + call ice_timer_start(timer_total) ! start timing entire run + call init_grid2 ! grid variables + call init_zbgc ! vertical biogeochemistry initialization + call init_calendar ! initialize some calendar stuff + call init_hist (dt) ! initialize output history file + + if (kdyn == 1) then + call init_evp + else if (kdyn == 2) then + call init_eap ! define eap dynamics parameters, variables + else if (kdyn == 3) then + call init_vp ! define vp dynamics parameters, variables + endif + + call init_coupler_flux ! initialize fluxes exchanged with coupler + + call init_thermo_vertical ! initialize vertical thermodynamics + + call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution + if (my_task == master_task) then + call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output + endif + + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution + floe_rad_l, & ! fsd size lower bound in m (radius) + floe_rad_c, & ! fsd size bin centre in m (radius) + floe_binwidth, & ! fsd size bin width in m (radius) + c_fsd_range, & ! string for history output + write_diags=(my_task == master_task)) ! write diag on master only + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call init_forcing_ocn(dt) ! initialize sss and sst from data + call init_state ! initialize the ice state + call init_transport ! initialize horizontal transport + call ice_HaloRestore_init ! restored boundary conditions + + call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & + wave_spec_out=wave_spec) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays + + call init_restart ! initialize restart variables + call init_diags ! initialize diagnostic output points + call init_history_therm ! initialize thermo history variables + call init_history_dyn ! initialize dynamic history variables + call calc_timesteps ! update timestep counter if not using npt_unit="1" + + call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) + call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_aero .or. tr_zaero) call faero_optics !initialize aerosol optical + !property tables + + ! Initialize shortwave components using swdn from previous timestep + ! if restarting. These components will be scaled to current forcing + ! in prep_radiation. + if (trim(runtype) == 'continue' .or. restart) & + call init_shortwave ! initialize radiative transfer + +! tcraig, use advance_timestep here +! istep = istep + 1 ! update time step counters +! istep1 = istep1 + 1 +! time = time + dt ! determine the time and date +! call calendar(time) ! at the end of the first timestep + call advance_timestep() + + !-------------------------------------------------------------------- + ! coupler communication or forcing data initialization + !-------------------------------------------------------------------- + + call init_forcing_atmo ! initialize atmospheric forcing (standalone) + + if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice + call get_forcing_atmo ! atmospheric forcing from data + call get_forcing_ocn(dt) ! ocean forcing from data + + ! isotopes + if (tr_iso) call fiso_default ! default values + ! aerosols + ! if (tr_aero) call faero_data ! data file + ! if (tr_zaero) call fzaero_data ! data file (gx1) + if (tr_aero .or. tr_zaero) call faero_default ! default values + if (skl_bgc .or. z_tracers) call get_forcing_bgc ! biogeochemistry + if (z_tracers) call get_atm_bgc ! biogeochemistry + + if (runtype == 'initial' .and. .not. restart) & + call init_shortwave ! initialize radiative transfer using current swdn + + call init_flux_atm ! initialize atmosphere fluxes sent to coupler + call init_flux_ocn ! initialize ocean fluxes sent to coupler + + if (write_ic) call accum_hist(dt) ! write initial conditions + + end subroutine cice_init + +!======================================================================= + + subroutine init_restart + + use ice_arrays_column, only: dhsn + use ice_blocks, only: nx_block, ny_block + use ice_calendar, only: calendar + use ice_constants, only: c0 + use ice_domain, only: nblocks + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_dyn_eap, only: read_restart_eap + use ice_dyn_shared, only: kdyn + use ice_grid, only: tmask + use ice_init, only: ice_ic + use ice_init_column, only: init_age, init_FY, init_lvl, & + init_meltponds_lvl, init_meltponds_topo, & + init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd + use ice_restart_column, only: restart_age, read_restart_age, & + restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & + restart_pond_lvl, read_restart_pond_lvl, & + restart_pond_topo, read_restart_pond_topo, & + restart_fsd, read_restart_fsd, & + restart_iso, read_restart_iso, & + restart_aero, read_restart_aero, & + restart_hbrine, read_restart_hbrine, & + restart_zsal, restart_bgc + use ice_restart_driver, only: restartfile + use ice_restart_shared, only: runtype, restart + use ice_state ! almost everything + + integer(kind=int_kind) :: & + i, j , & ! horizontal indices + iblk ! block index + logical(kind=log_kind) :: & + tr_iage, tr_FY, tr_lvl, tr_pond_lvl, & + tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + skl_bgc, z_tracers, solve_zsal + integer(kind=int_kind) :: & + ntrcr + integer(kind=int_kind) :: & + nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice + + character(len=*), parameter :: subname = '(init_restart)' + + call icepack_query_tracer_sizes(ntrcr_out=ntrcr) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call icepack_query_parameters(skl_bgc_out=skl_bgc, & + z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_lvl_out=tr_lvl, tr_pond_lvl_out=tr_pond_lvl, & + tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (trim(runtype) == 'continue') then + ! start from core restart file + call restartfile() ! given by pointer in ice_in + call calendar() ! update time parameters + if (kdyn == 2) call read_restart_eap ! EAP + else if (restart) then ! ice_ic = core restart file + call restartfile (ice_ic) ! or 'default' or 'none' + !!! uncomment to create netcdf + ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file + !!! uncomment if EAP restart data exists + ! if (kdyn == 2) call read_restart_eap + endif + + ! tracers + ! ice age tracer + if (tr_iage) then + if (trim(runtype) == 'continue') & + restart_age = .true. + if (restart_age) then + call read_restart_age + else + do iblk = 1, nblocks + call init_age(trcrn(:,:,nt_iage,:,iblk)) + enddo ! iblk + endif + endif + ! first-year area tracer + if (tr_FY) then + if (trim(runtype) == 'continue') restart_FY = .true. + if (restart_FY) then + call read_restart_FY + else + do iblk = 1, nblocks + call init_FY(trcrn(:,:,nt_FY,:,iblk)) + enddo ! iblk + endif + endif + ! level ice tracer + if (tr_lvl) then + if (trim(runtype) == 'continue') restart_lvl = .true. + if (restart_lvl) then + call read_restart_lvl + else + do iblk = 1, nblocks + call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & + trcrn(:,:,nt_vlvl,:,iblk)) + enddo ! iblk + endif + endif + ! level-ice melt ponds + if (tr_pond_lvl) then + if (trim(runtype) == 'continue') & + restart_pond_lvl = .true. + if (restart_pond_lvl) then + call read_restart_pond_lvl + else + do iblk = 1, nblocks + call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk), & + dhsn(:,:,:,iblk)) + enddo ! iblk + endif + endif + ! topographic melt ponds + if (tr_pond_topo) then + if (trim(runtype) == 'continue') & + restart_pond_topo = .true. + if (restart_pond_topo) then + call read_restart_pond_topo + else + do iblk = 1, nblocks + call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk)) + enddo ! iblk + endif ! .not. restart_pond + endif + ! floe size distribution + if (tr_fsd) then + if (trim(runtype) == 'continue') restart_fsd = .true. + if (restart_fsd) then + call read_restart_fsd + else + call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) + endif + endif + + ! isotopes + if (tr_iso) then + if (trim(runtype) == 'continue') restart_iso = .true. + if (restart_iso) then + call read_restart_iso + else + do iblk = 1, nblocks + call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & + trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) + enddo ! iblk + endif + endif + + if (tr_aero) then ! ice aerosol + if (trim(runtype) == 'continue') restart_aero = .true. + if (restart_aero) then + call read_restart_aero + else + do iblk = 1, nblocks + call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) + enddo ! iblk + endif ! .not. restart_aero + endif + + if (trim(runtype) == 'continue') then + if (tr_brine) & + restart_hbrine = .true. + if (solve_zsal) & + restart_zsal = .true. + if (skl_bgc .or. z_tracers) & + restart_bgc = .true. + endif + + if (tr_brine .or. skl_bgc) then ! brine height tracer + call init_hbrine + if (tr_brine .and. restart_hbrine) call read_restart_hbrine + endif + + if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry + if (tr_fsd) then + write (nu_diag,*) 'FSD implementation incomplete for use with BGC' + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + endif + call init_bgc + endif + + !----------------------------------------------------------------- + ! aggregate tracers + !----------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + if (tmask(i,j,iblk)) then + call icepack_aggregate(ncat = ncat, & + aicen = aicen(i,j,:,iblk), & + trcrn = trcrn(i,j,:,:,iblk), & + vicen = vicen(i,j,:,iblk), & + vsnon = vsnon(i,j,:,iblk), & + aice = aice (i,j, iblk), & + trcr = trcr (i,j,:,iblk), & + vice = vice (i,j, iblk), & + vsno = vsno (i,j, iblk), & + aice0 = aice0(i,j, iblk), & + ntrcr = ntrcr, & + trcr_depend = trcr_depend, & + trcr_base = trcr_base, & + n_trcr_strata = n_trcr_strata, & + nt_strata = nt_strata) + else + ! tcraig, reset all tracer values on land to zero + trcrn(i,j,:,:,iblk) = c0 + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + end subroutine init_restart + +!======================================================================= + + end module CICE_InitMod + +!======================================================================= diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 new file mode 100644 index 000000000..231fc9dd3 --- /dev/null +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -0,0 +1,608 @@ + + module halochk_data + + use CICE_InitMod + use ice_kinds_mod, only: int_kind, dbl_kind, real_kind, log_kind + use ice_blocks, only: block, get_block, nx_block, ny_block, nblocks_tot, nghost + use ice_boundary, only: ice_HaloUpdate + use ice_constants, only: c0, c1, & + field_loc_unknown, field_loc_noupdate, & + field_loc_center, field_loc_NEcorner, & + field_loc_Nface, field_loc_Eface, & + field_type_unknown, field_type_noupdate, & + field_type_scalar, field_type_vector, field_type_angle + use ice_communicate, only: my_task, master_task, get_num_procs, MPI_COMM_ICE + use ice_distribution, only: ice_distributionGetBlockID, ice_distributionGet + use ice_domain_size, only: nx_global, ny_global, & + block_size_x, block_size_y, max_blocks + use ice_domain, only: distrb_info, halo_info, & + ew_boundary_type, ns_boundary_type + use ice_exit, only: abort_ice, end_run + use ice_global_reductions, only: global_minval, global_maxval, global_sum + + implicit none + + integer(int_kind), parameter :: & + passflag = 0, & + failflag = 1 + + end module halochk_data + + + program halochk + + ! This tests the CICE halo update methods by + ! using CICE_InitMod (from the standalone model) to read/initialize + ! a CICE grid/configuration. + + use halochk_data + + implicit none + + integer(int_kind) :: nn, nl, nt, i, j, k1, k2, n, ib, ie, jb, je, iblock + integer(int_kind) :: blockID, numBlocks + type (block) :: this_block + + real(dbl_kind) , allocatable :: darrayi1(:,:,:) , darrayj1(:,:,:) + real(dbl_kind) , allocatable :: darrayi2(:,:,:,:) , darrayj2(:,:,:,:) + real(dbl_kind) , allocatable :: darrayi3(:,:,:,:,:), darrayj3(:,:,:,:,:) + real(real_kind) , allocatable :: rarrayi1(:,:,:) , rarrayj1(:,:,:) + real(real_kind) , allocatable :: rarrayi2(:,:,:,:) , rarrayj2(:,:,:,:) + real(real_kind) , allocatable :: rarrayi3(:,:,:,:,:), rarrayj3(:,:,:,:,:) + integer(int_kind), allocatable :: iarrayi1(:,:,:) , iarrayj1(:,:,:) + integer(int_kind), allocatable :: iarrayi2(:,:,:,:) , iarrayj2(:,:,:,:) + integer(int_kind), allocatable :: iarrayi3(:,:,:,:,:), iarrayj3(:,:,:,:,:) + logical(log_kind), allocatable :: larrayi1(:,:,:) , larrayj1(:,:,:) + + real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) + real(dbl_kind), allocatable :: cidata_nup(:,:,:,:,:),cjdata_nup(:,:,:,:,:) + real(dbl_kind), allocatable :: cidata_std(:,:,:,:,:),cjdata_std(:,:,:,:,:) + + integer(int_kind), parameter :: maxtests = 10 + integer(int_kind), parameter :: maxtypes = 4 + integer(int_kind), parameter :: maxlocs = 5 + integer(int_kind), parameter :: nz1 = 3 + integer(int_kind), parameter :: nz2 = 4 + real(dbl_kind) :: aichk,ajchk,cichk,cjchk + character(len=16) :: locs_name(maxlocs), types_name(maxtypes) + integer(int_kind) :: field_loc(maxlocs), field_type(maxtypes) + integer(int_kind) :: npes, ierr, ntask, testcnt, tottest + integer(int_kind) :: errorflag0, gflag, k1m, k2m, ptcntsum + integer(int_kind), allocatable :: errorflag(:) + integer(int_kind), allocatable :: ptcnt(:) + character(len=128), allocatable :: stringflag(:) + character(len=32) :: string + logical :: first_call = .true. + + real(dbl_kind) , parameter :: fillval = -88888.0_dbl_kind + real(dbl_kind) , parameter :: dhalofillval = -999.0_dbl_kind + real(real_kind) , parameter :: rhalofillval = -999.0_real_kind + integer(int_kind), parameter :: ihalofillval = -999 + character(len=*) , parameter :: subname='(halochk)' + + !----------------------------------------------------------------- + ! Initialize CICE + !----------------------------------------------------------------- + + call CICE_Initialize + npes = get_num_procs() + + locs_name (:) = 'unknown' + types_name(:) = 'unknown' + field_type(:) = field_type_unknown + field_loc (:) = field_loc_unknown + + types_name(1) = 'scalar' + field_type(1) = field_type_scalar + types_name(2) = 'vector' + field_type(2) = field_type_vector + types_name(3) = 'angle' + field_type(3) = field_type_angle + types_name(4) = 'none' + field_type(4) = field_type_noupdate +! types_name(5) = 'unknown' +! field_type(5) = field_type_unknown ! aborts in CICE, as expected + + locs_name (1) = 'center' + field_loc (1) = field_loc_center + locs_name (2) = 'NEcorner' + field_loc (2) = field_loc_NEcorner + locs_name (3) = 'Nface' + field_loc (3) = field_loc_Nface + locs_name (4) = 'Eface' + field_loc (4) = field_loc_Eface + locs_name (5) = 'none' + field_loc (5) = field_loc_noupdate +! locs_name (6) = 'unknown' +! field_loc (6) = field_loc_unknown ! aborts in CICE, as expected + + tottest = maxtests * maxlocs * maxtypes + allocate(errorflag(tottest)) + allocate(stringflag(tottest)) + allocate(ptcnt(tottest)) + ptcnt(:) = 0 + + !----------------------------------------------------------------- + ! Testing + !----------------------------------------------------------------- + + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) '==========================================================' + write(6,*) ' ' + write(6,*) 'RunningUnitTest HALOCHK' + write(6,*) ' ' + write(6,*) ' npes = ',npes + write(6,*) ' my_task = ',my_task + write(6,*) ' nx_global = ',nx_global + write(6,*) ' ny_global = ',ny_global + write(6,*) ' block_size_x = ',block_size_x + write(6,*) ' block_size_y = ',block_size_y + write(6,*) ' nblocks_tot = ',nblocks_tot + write(6,*) ' tottest = ',tottest + write(6,*) ' ' + endif + + errorflag0 = passflag + errorflag(:) = passflag + stringflag(:) = ' ' + + ! --------------------------- + ! TEST HALO UPDATE + ! --------------------------- + + if (my_task == master_task) write(6,*) ' ' + + allocate(darrayi1 (nx_block,ny_block,max_blocks)) + allocate(darrayj1 (nx_block,ny_block,max_blocks)) + allocate(darrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(darrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(darrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(darrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(rarrayi1 (nx_block,ny_block,max_blocks)) + allocate(rarrayj1 (nx_block,ny_block,max_blocks)) + allocate(rarrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(rarrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(rarrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(rarrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(iarrayi1 (nx_block,ny_block,max_blocks)) + allocate(iarrayj1 (nx_block,ny_block,max_blocks)) + allocate(iarrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(iarrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(iarrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(iarrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(larrayi1 (nx_block,ny_block,max_blocks)) + allocate(larrayj1 (nx_block,ny_block,max_blocks)) + allocate(cidata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cidata_std(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_std(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cidata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) + + darrayi1 = fillval + darrayj1 = fillval + darrayi2 = fillval + darrayj2 = fillval + darrayi3 = fillval + darrayj3 = fillval + rarrayi1 = fillval + rarrayj1 = fillval + rarrayi2 = fillval + rarrayj2 = fillval + rarrayi3 = fillval + rarrayj3 = fillval + iarrayi1 = fillval + iarrayj1 = fillval + iarrayi2 = fillval + iarrayj2 = fillval + iarrayi3 = fillval + iarrayj3 = fillval + larrayi1 = .false. + larrayj1 = .true. + cidata_bas = fillval + cjdata_bas = fillval + cidata_std = fillval + cjdata_std = fillval + cidata_nup = fillval + cjdata_nup = fillval + + call ice_distributionGet(distrb_info, numLocalBlocks = numBlocks) + + !--- baseline data --- + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + do k2 = 1,nz2 + do k1 = 1,nz1 + do j = 1,ny_block + do i = 1,nx_block + cidata_bas(i,j,k1,k2,iblock) = real(this_block%i_glob(i),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind + cjdata_bas(i,j,k1,k2,iblock) = real(this_block%j_glob(j),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind + enddo + enddo + enddo + enddo + enddo + + !--- setup nup (noupdate) solution, set halo/pad will fillval --- + + cidata_nup(:,:,:,:,:) = cidata_bas(:,:,:,:,:) + cjdata_nup(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + cidata_nup(1:ib-1 ,: ,:,:,iblock) = fillval + cjdata_nup(1:ib-1 ,: ,:,:,iblock) = fillval + cidata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval + cjdata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval + cidata_nup(: ,1:jb-1 ,:,:,iblock) = fillval + cjdata_nup(: ,1:jb-1 ,:,:,iblock) = fillval + cidata_nup(: ,je+1:ny_block,:,:,iblock) = fillval + cjdata_nup(: ,je+1:ny_block,:,:,iblock) = fillval + enddo + + !--- setup std solution for cyclic, closed, open, tripole solution --- + + cidata_std(:,:,:,:,:) = cidata_bas(:,:,:,:,:) + cjdata_std(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) + + !--- halo off on east and west boundary --- + if (ew_boundary_type == 'closed' .or. & + ew_boundary_type == 'open' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%i_glob(ib) == 1) then + cidata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval + cjdata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval + endif + if (this_block%i_glob(ie) == nx_global) then + cidata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval + cjdata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--- halo off on south boundary --- + if (ns_boundary_type == 'closed' .or. & + ns_boundary_type == 'open' .or. & + ns_boundary_type == 'tripole' .or. & + ns_boundary_type == 'tripoleT' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%j_glob(jb) == 1) then + cidata_std(:,1:jb-1,:,:,iblock) = dhalofillval + cjdata_std(:,1:jb-1,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--- halo off on north boundary --- + if (ns_boundary_type == 'closed' .or. & + ns_boundary_type == 'open' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%j_glob(je) == ny_global) then + cidata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval + cjdata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--------------------------------------------------------------- + + testcnt = 0 + do nn = 1, maxtests + do nl = 1, maxlocs + do nt = 1, maxtypes + + !--- setup test --- + first_call = .true. + testcnt = testcnt + 1 + if (testcnt > tottest) then + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK FAILED' + write(6,*) ' ' + call abort_ice(subname//' testcnt > tottest',file=__FILE__,line=__LINE__) + endif + endif + + !--- fill arrays --- + darrayi1(:,:,:) = fillval + darrayj1(:,:,:) = fillval + darrayi2(:,:,:,:) = fillval + darrayj2(:,:,:,:) = fillval + darrayi3(:,:,:,:,:) = fillval + darrayj3(:,:,:,:,:) = fillval + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + do j = jb,je + do i = ib,ie + k1 = 1 + k2 = 1 + darrayi1(i,j,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj1(i,j,iblock) = cjdata_bas(i,j,k1,k2,iblock) + do k1 = 1,nz1 + k2 = 1 + darrayi2(i,j,k1,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj2(i,j,k1,iblock) = cjdata_bas(i,j,k1,k2,iblock) + do k2 = 1,nz2 + darrayi3(i,j,k1,k2,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj3(i,j,k1,k2,iblock) = cjdata_bas(i,j,k1,k2,iblock) + enddo + enddo + enddo + enddo + enddo + + + !--- halo update --- + + if (nn == 1) then + k1m = 1 + k2m = 1 + string = '2DR8' + call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 2) then + k1m = nz1 + k2m = 1 + string = '3DR8' + call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 3) then + k1m = nz1 + k2m = nz2 + string = '4DR8' + call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 4) then + k1m = 1 + k2m = 1 + string = '2DR4' + rarrayi1 = real(darrayi1,kind=real_kind) + rarrayj1 = real(darrayj1,kind=real_kind) + call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi1 = real(rarrayi1,kind=dbl_kind) + darrayj1 = real(rarrayj1,kind=dbl_kind) + elseif (nn == 5) then + k1m = nz1 + k2m = 1 + string = '3DR4' + rarrayi2 = real(darrayi2,kind=real_kind) + rarrayj2 = real(darrayj2,kind=real_kind) + call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi2 = real(rarrayi2,kind=dbl_kind) + darrayj2 = real(rarrayj2,kind=dbl_kind) + elseif (nn == 6) then + k1m = nz1 + k2m = nz2 + string = '4DR4' + rarrayi3 = real(darrayi3,kind=real_kind) + rarrayj3 = real(darrayj3,kind=real_kind) + call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi3 = real(rarrayi3,kind=dbl_kind) + darrayj3 = real(rarrayj3,kind=dbl_kind) + elseif (nn == 7) then + k1m = 1 + k2m = 1 + string = '2DI4' + iarrayi1 = nint(darrayi1) + iarrayj1 = nint(darrayj1) + call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi1 = real(iarrayi1,kind=dbl_kind) + darrayj1 = real(iarrayj1,kind=dbl_kind) + elseif (nn == 8) then + k1m = nz1 + k2m = 1 + string = '3DI4' + iarrayi2 = nint(darrayi2) + iarrayj2 = nint(darrayj2) + call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi2 = real(iarrayi2,kind=dbl_kind) + darrayj2 = real(iarrayj2,kind=dbl_kind) + elseif (nn == 9) then + k1m = nz1 + k2m = nz2 + string = '4DI4' + iarrayi3 = nint(darrayi3) + iarrayj3 = nint(darrayj3) + call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi3 = real(iarrayi3,kind=dbl_kind) + darrayj3 = real(iarrayj3,kind=dbl_kind) + elseif (nn == 10) then + k1m = 1 + k2m = 1 + string = '2DL1' + larrayi1 = .true. + where (darrayi1 == fillval) larrayi1 = .false. + larrayj1 = .false. + where (darrayj1 == fillval) larrayj1 = .true. + call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=0) + call ice_haloUpdate(larrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=1) + darrayi1 = c0 + where (larrayi1) darrayi1 = c1 + darrayj1 = c0 + where (larrayj1) darrayj1 = c1 + endif + + write(stringflag(testcnt),'(5a10)') trim(string),trim(locs_name(nl)),trim(types_name(nt)), & + trim(ew_boundary_type),trim(ns_boundary_type) + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + ! just check non-padded gridcells +! do j = 1,ny_block +! do i = 1,nx_block + do j = jb-nghost, je+nghost + do i = ib-nghost, ie+nghost + do k1 = 1,k1m + do k2 = 1,k2m + if (nn == 1 .or. nn == 4 .or. nn == 7 .or. nn == 10) then + aichk = darrayi1(i,j,iblock) + ajchk = darrayj1(i,j,iblock) + elseif (nn == 2 .or. nn == 5 .or. nn == 8) then + aichk = darrayi2(i,j,k1,iblock) + ajchk = darrayj2(i,j,k1,iblock) + elseif (nn == 3 .or. nn == 6 .or. nn == 9) then + aichk = darrayi3(i,j,k1,k2,iblock) + ajchk = darrayj3(i,j,k1,k2,iblock) + endif + + if (field_loc (nl) == field_loc_noupdate .or. & + field_type(nt) == field_type_noupdate) then + cichk = cidata_nup(i,j,k1,k2,iblock) + cjchk = cjdata_nup(i,j,k1,k2,iblock) + else + cichk = cidata_std(i,j,k1,k2,iblock) + cjchk = cjdata_std(i,j,k1,k2,iblock) + endif + + if (nn >= 7 .and. nn <= 9) then + cichk = real(nint(cichk),kind=dbl_kind) + cjchk = real(nint(cjchk),kind=dbl_kind) + endif + + if (nn == 10) then + if (cichk == dhalofillval .or. cichk == fillval) then + cichk = c0 + else + cichk = c1 + endif + if (cjchk == dhalofillval .or. cjchk == fillval) then + cjchk = c1 + else + cjchk = c0 + endif + endif + + ptcnt(testcnt) = ptcnt(testcnt) + 1 + call chkresults(aichk,cichk,errorflag(testcnt),testcnt, & + i,j,k1,k2,iblock,first_call,stringflag(testcnt),trim(string)//'_I') + call chkresults(ajchk,cjchk,errorflag(testcnt),testcnt, & + i,j,k1,k2,iblock,first_call,stringflag(testcnt),trim(string)//'_J') + enddo ! k2 + enddo ! k1 + enddo ! i + enddo ! j + enddo ! iblock + + enddo ! maxtypes + enddo ! maxlocs + enddo ! maxtests + + ! --------------------------- + ! SUMMARY + ! --------------------------- + + do n = 1,tottest + gflag = global_maxval(errorflag(n), MPI_COMM_ICE) + errorflag(n) = gflag + ptcntsum = global_sum(ptcnt(n),distrb_info) + ptcnt(n) = ptcntsum + enddo + errorflag0 = maxval(errorflag(:)) + + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK COMPLETED SUCCESSFULLY' + write(6,*) ' ' + do n = 1,tottest + if (errorflag(n) == passflag) then + write(6,*) 'PASS ',trim(stringflag(n)),ptcnt(n) + else + write(6,*) 'FAIL ',trim(stringflag(n)),ptcnt(n) + endif + enddo + write(6,*) ' ' + if (errorflag0 == passflag) then + write(6,*) 'HALOCHK TEST COMPLETED SUCCESSFULLY' + else + write(6,*) 'HALOCHK TEST FAILED' + endif + write(6,*) ' ' + write(6,*) '==========================================================' + write(6,*) ' ' + endif + + + !----------------------------------------------------------------- + ! Gracefully end + !----------------------------------------------------------------- + + call end_run() + + end program halochk + +!======================================================================= + + subroutine chkresults(a1,r1,errorflag,testcnt,i,j,k1,k2,iblock,first_call,stringflag,string) + + use halochk_data + + implicit none + + real(dbl_kind) , intent(in) :: a1,r1 + integer(int_kind), intent(inout) :: errorflag + integer(int_kind), intent(in) :: i,j,k1,k2,iblock,testcnt + logical , intent(inout) :: first_call + character(len=*) , intent(in) :: stringflag,string + + logical,parameter :: print_always = .false. + character(len=*) , parameter :: subname='(chkresults)' + + if (a1 /= r1 .or. print_always) then + errorflag = failflag + if (first_call) then + write(100+my_task,*) ' ' + write(100+my_task,'(a,i4,2a)') '------- TEST = ',testcnt,' ',trim(stringflag) + write(100+my_task,'(a)') ' test task i j k1 k2 iblock expected computed diff' + first_call = .false. + endif + write(100+my_task,1001) trim(string),testcnt,my_task,i,j,k1,k2,iblock,r1,a1,r1-a1 + endif + + 1001 format(a8,7i6,3f12.3) + + end subroutine chkresults +!======================================================================= diff --git a/cicecore/shared/ice_constants.F90 b/cicecore/shared/ice_constants.F90 index f2da2ef9d..6656213be 100644 --- a/cicecore/shared/ice_constants.F90 +++ b/cicecore/shared/ice_constants.F90 @@ -97,8 +97,7 @@ module ice_constants field_loc_center = 1, & field_loc_NEcorner = 2, & field_loc_Nface = 3, & - field_loc_Eface = 4, & - field_loc_Wface = 5 + field_loc_Eface = 4 !----------------------------------------------------------------- ! field type attribute - necessary for handling diff --git a/configuration/scripts/Makefile b/configuration/scripts/Makefile index a2f17256f..872f426ad 100644 --- a/configuration/scripts/Makefile +++ b/configuration/scripts/Makefile @@ -74,7 +74,7 @@ AR := ar .SUFFIXES: -.PHONY: all cice libcice targets target db_files db_flags clean realclean helloworld calchk sumchk bcstchk gridavgchk optargs +.PHONY: all cice libcice targets target db_files db_flags clean realclean helloworld calchk sumchk bcstchk gridavgchk halochk optargs all: $(EXEC) cice: $(EXEC) @@ -93,7 +93,7 @@ targets: @echo " " @echo "Supported Makefile Targets are: cice, libcice, makdep, depends, clean, realclean" @echo " Diagnostics: targets, db_files, db_flags" - @echo " Unit Tests : helloworld, calchk, sumchk, bcstchk, gridavgchk, optargs" + @echo " Unit Tests : helloworld, calchk, sumchk, bcstchk, gridavgchk, halochk, optargs" target: targets db_files: @@ -151,6 +151,8 @@ bcstchk: $(EXEC) gridavgchk: $(EXEC) +halochk: $(EXEC) + # this builds just a subset of source code specified explicitly and requires a separate target HWOBJS := helloworld.o diff --git a/configuration/scripts/cice_decomp.csh b/configuration/scripts/cice_decomp.csh index 0c6715f3b..bcf27beee 100755 --- a/configuration/scripts/cice_decomp.csh +++ b/configuration/scripts/cice_decomp.csh @@ -156,8 +156,10 @@ if (${ICE_DECOMP_MXBLCKS} > 0) set mxblcks = ${ICE_DECOMP_MXBLCKS} set decomp = 'cartesian' set dshape = 'slenderX2' -if (${nxglob} % ${cicepes} != 0) set decomp = 'roundrobin' -if (${mxblcks} * ${blcky} * 2 < ${nyglob}) set decomp = 'roundrobin' +if (${cicepes} % 2 != 0) set decomp = 'roundrobin' +if (${nyglob} % (${blcky} * 2) != 0) set decomp = 'roundrobin' +if (${nxglob} % ${blckx} != 0) set decomp = 'roundrobin' +if (((${nxglob} * 2) % (${cicepes} * ${blckx})) != 0) set decomp = 'roundrobin' #--- outputs --- diff --git a/configuration/scripts/options/set_env.halochk b/configuration/scripts/options/set_env.halochk new file mode 100644 index 000000000..d09166d2f --- /dev/null +++ b/configuration/scripts/options/set_env.halochk @@ -0,0 +1,2 @@ +setenv ICE_DRVOPT unittest/halochk +setenv ICE_TARGET halochk diff --git a/configuration/scripts/options/set_nml.bcyclic b/configuration/scripts/options/set_nml.bcyclic new file mode 100644 index 000000000..3a5ae1a7b --- /dev/null +++ b/configuration/scripts/options/set_nml.bcyclic @@ -0,0 +1,3 @@ +ew_boundary_type = 'cyclic' +ns_boundary_type = 'cyclic' + diff --git a/configuration/scripts/options/set_nml.bopen b/configuration/scripts/options/set_nml.bopen new file mode 100644 index 000000000..0e2d5f388 --- /dev/null +++ b/configuration/scripts/options/set_nml.bopen @@ -0,0 +1,3 @@ +ew_boundary_type = 'open' +ns_boundary_type = 'open' + diff --git a/doc/source/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index 289f626a9..606ae1397 100644 --- a/doc/source/user_guide/ug_testing.rst +++ b/doc/source/user_guide/ug_testing.rst @@ -736,6 +736,7 @@ The following are brief descriptions of some of the current unit tests, - **calchk** is a unit test that exercises the CICE calendar over 100,000 years and verifies correctness. This test does not depend on the CICE initialization. - **gridavgchk** is a unit test that exercises the CICE grid_average_X2Y methods and verifies results. + - **halochk** is a unit test that exercises the CICE haloUpdate methods and verifies results. - **helloworld** is a simple test that writes out helloworld and uses no CICE infrastructure. This tests exists to demonstrate how to build a unit test by specifying the object files directly in the Makefile From c1148c3448939b2fb879517fd122e3f234c99c4e Mon Sep 17 00:00:00 2001 From: apcraig Date: Thu, 9 Mar 2023 18:28:00 -0700 Subject: [PATCH 02/13] - Fix bug in serial/ice_boundary.F90 tripoleT halo update - Reduce redundant tripole buffer copies in serial/ice_boundary.F90 - Generalize iSrc wraparound calculation in ice_boundary.F90 - Add open, cyclic, tripole, and tripoleT set_nml files - Update unittest suite --- .../infrastructure/comm/mpi/ice_boundary.F90 | 58 ++++- .../comm/serial/ice_boundary.F90 | 185 +++++++------- cicecore/drivers/unittest/halochk/halochk.F90 | 241 ++++++++++++++---- .../{set_nml.bcyclic => set_nml.cyclic} | 0 .../options/{set_nml.bopen => set_nml.open} | 0 configuration/scripts/options/set_nml.tripole | 3 + .../scripts/options/set_nml.tripolet | 3 + configuration/scripts/tests/unittest_suite.ts | 41 ++- 8 files changed, 379 insertions(+), 152 deletions(-) rename configuration/scripts/options/{set_nml.bcyclic => set_nml.cyclic} (100%) rename configuration/scripts/options/{set_nml.bopen => set_nml.open} (100%) create mode 100644 configuration/scripts/options/set_nml.tripole create mode 100644 configuration/scripts/options/set_nml.tripolet diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 index 77af8ceee..58396a58e 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 @@ -12,6 +12,43 @@ module ice_boundary ! fixes for non-existent blocks ! 2008-01-28: Elizabeth Hunke replaced old routines with new POP ! infrastructure +! +!----------------------------------------------------------------------- +! +! Some notes on tripole, A-H below are gridpoints at i = 1:nx_global +! where nx_global=8. The schematics below show the general layout of the center +! points on the tripole fold. More complex pictures are needed to show +! relative orientation and offsets of east, north, and northeast points +! across the fold. See also appendix E of the NEMO_manual, +! https://zenodo.org/record/6334656#.YiYirhPMLXQ. Note the NFtype=T +! is the tripole u-fold grid with T-grid=center, U-grid=east, V-grid=north, +! and F-grid=northeast points in CICE. NFtype=F is similar to tripoleT +! except for the treatment of the poles. The CICE implementation also +! averages all degenerate points, NEMO's strategy seems to be to copy +! data from one side of the tripole to the other for degenerate points. +! +! tripole: u-fold, fold is on north edge of ny_global +! north and northeast points on the fold are degenerate and averaged +! A,H,D,and E are pole points +! +! ny_global+2 H G F E D C B A @ny_global-1 +! ny_global+1 H G F E D C B A @ny_global +! ny_global A B C D E F G H +! ny_global-1 A B C D E F G H +! +! tripoleT: t-fold, fold is thru center of ny_global +! center and east points at ny_global are degenerate and averaged +! north and northeast point at ny_global are not prognostic, they are halos +! A and E are pole points +! +! ny_global+2 H G F E D C B A @ny_global-2 +! ny_global+1 H G F E D C B A @ny_global-1 +! ny_global A BH CG DF E FD GC HB A +! ny_global-1 A B C D E F G H +! ny_global-2 A B C D E F G H +! +!----------------------------------------------------------------------- + use mpi ! MPI Fortran module use ice_kinds_mod @@ -1571,7 +1608,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1986,7 +2023,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2401,7 +2438,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2945,7 +2982,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3419,7 +3456,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3893,7 +3930,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -4389,7 +4426,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -4887,7 +4924,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -5385,7 +5422,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -5720,7 +5757,8 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal + if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface, do not need to replace !*** top row of physical domain, so jSrc should be diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index 65da2a6f4..30a513430 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -12,6 +12,11 @@ module ice_boundary ! fixes for non-existent blocks ! 2008-01-28: Elizabeth Hunke replaced old routines with new POP ! infrastructure +! 2023-03-09: Tony Craig updated the implementation to fix bug in +! tripoleT and reduce number of copies in tripole overall. +! Because all blocks are local, can fill the tripole +! buffer from "north" copies. This is not true for +! the MPI version. use ice_kinds_mod use ice_communicate, only: my_task @@ -312,10 +317,11 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** for that !echmod if (tripoleBlock .and. dstProc /= srcProc) then - if (tripoleBlock) then - call ice_HaloIncrementMsgCount(sendCount, recvCount, & - srcProc, dstProc, northMsgSize) - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloIncrementMsgCount(sendCount, recvCount, & +! srcProc, dstProc, northMsgSize) +! endif !*** find west neighbor block and add to message count @@ -338,10 +344,11 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** for that !echmod if (tripoleBlock .and. dstProc /= srcProc) then - if (tripoleBlock) then - call ice_HaloIncrementMsgCount(sendCount, recvCount, & - srcProc, dstProc, northMsgSize) - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloIncrementMsgCount(sendCount, recvCount, & +! srcProc, dstProc, northMsgSize) +! endif !*** find northeast neighbor block and add to message count @@ -354,11 +361,12 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_distributionGetBlockLoc(dist, neBlock, dstProc, & dstLocalID) - else if (neBlock < 0) then ! tripole north row - msgSize = northMsgSize ! tripole needs whole top row of block - - call ice_distributionGetBlockLoc(dist, abs(neBlock), dstProc, & - dstLocalID) +! tcx,tcraig, 3/2023, this is not needed +! else if (neBlock < 0) then ! tripole north row +! msgSize = northMsgSize ! tripole needs whole top row of block +! +! call ice_distributionGetBlockLoc(dist, abs(neBlock), dstProc, & +! dstLocalID) else dstProc = 0 dstLocalID = 0 @@ -378,11 +386,12 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_distributionGetBlockLoc(dist, nwBlock, dstProc, & dstLocalID) - else if (nwBlock < 0) then ! tripole north row, count block - msgSize = northMsgSize ! tripole NE corner update - entire row needed - - call ice_distributionGetBlockLoc(dist, abs(nwBlock), dstProc, & - dstLocalID) +! tcx,tcraig, 3/2023, this is not needed +! else if (nwBlock < 0) then ! tripole north row, count block +! msgSize = northMsgSize ! tripole NE corner update - entire row needed +! +! call ice_distributionGetBlockLoc(dist, abs(nwBlock), dstProc, & +! dstLocalID) else dstProc = 0 @@ -484,8 +493,6 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & northBlock = ice_blocksGetNbrID(iblock, ice_blocksNorth, & ewBoundaryType, nsBoundaryType) - call ice_HaloMsgCreate(halo, dist, iblock, northBlock, 'north') - !*** set tripole flag and add two copies for inserting !*** and extracting info from the tripole buffer @@ -495,6 +502,7 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_HaloMsgCreate(halo, dist, -iblock, iblock, 'north') else tripoleBlock = .false. + call ice_HaloMsgCreate(halo, dist, iblock, northBlock, 'north') endif !*** find south neighbor block @@ -515,9 +523,10 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** the east block to make sure enough information is !*** available for tripole manipulations - if (tripoleBlock) then - call ice_HaloMsgCreate(halo, dist, iblock, -eastBlock, 'north') - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloMsgCreate(halo, dist, iblock, -eastBlock, 'north') +! endif !*** find west neighbor block @@ -530,9 +539,10 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** the west block to make sure enough information is !*** available for tripole manipulations - if (tripoleBlock) then - call ice_HaloMsgCreate(halo, dist, iblock, -westBlock, 'north') - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloMsgCreate(halo, dist, iblock, -westBlock, 'north') +! endif !*** find northeast neighbor block @@ -955,7 +965,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1264,7 +1274,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1573,7 +1583,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1983,7 +1993,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2318,7 +2328,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2653,7 +2663,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3003,7 +3013,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3355,7 +3365,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3707,7 +3717,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3932,7 +3942,8 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal + if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface, do not need to replace !*** top row of physical domain, so jSrc should be @@ -4348,36 +4359,37 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) halo%numLocalCopies = msgIndx - else - - !*** tripole grid - copy entire top halo+1 - !*** rows into global buffer at src location - - msgIndx = halo%numLocalCopies - - do j=1,nghost+1 - do i=1,ieSrc-ibSrc+1 - - msgIndx = msgIndx + 1 - - halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 - halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j - halo%srcLocalAddr(3,msgIndx) = srcLocalID - - halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) - halo%dstLocalAddr(2,msgIndx) = j - halo%dstLocalAddr(3,msgIndx) = -dstLocalID - - end do - end do - - halo%numLocalCopies = msgIndx +! tcx,tcraig, 3/2023, this is not needed +! else +! +! !*** tripole grid - copy entire top halo+1 +! !*** rows into global buffer at src location +! +! msgIndx = halo%numLocalCopies +! +! do j=1,nghost+1 +! do i=1,ieSrc-ibSrc+1 +! +! msgIndx = msgIndx + 1 +! +! halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 +! halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j +! halo%srcLocalAddr(3,msgIndx) = srcLocalID +! +! halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) +! halo%dstLocalAddr(2,msgIndx) = j +! halo%dstLocalAddr(3,msgIndx) = -dstLocalID +! +! end do +! end do +! +! halo%numLocalCopies = msgIndx endif case ('northwest') - !*** normal northeast boundary - just copy NW corner + !*** normal northwest boundary - just copy NW corner !*** of physical domain into SE halo of NW nbr block if (dstBlock > 0) then @@ -4402,30 +4414,31 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) halo%numLocalCopies = msgIndx - else - - !*** tripole grid - copy entire top halo+1 - !*** rows into global buffer at src location - - msgIndx = halo%numLocalCopies - - do j=1,nghost+1 - do i=1,ieSrc-ibSrc+1 - - msgIndx = msgIndx + 1 - - halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 - halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j - halo%srcLocalAddr(3,msgIndx) = srcLocalID - - halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) - halo%dstLocalAddr(2,msgIndx) = j - halo%dstLocalAddr(3,msgIndx) = -dstLocalID - - end do - end do - - halo%numLocalCopies = msgIndx +! tcx,tcraig, 3/2023, this is not needed +! else +! +! !*** tripole grid - copy entire top halo+1 +! !*** rows into global buffer at src location +! +! msgIndx = halo%numLocalCopies +! +! do j=1,nghost+1 +! do i=1,ieSrc-ibSrc+1 +! +! msgIndx = msgIndx + 1 +! +! halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 +! halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j +! halo%srcLocalAddr(3,msgIndx) = srcLocalID +! +! halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) +! halo%dstLocalAddr(2,msgIndx) = j +! halo%dstLocalAddr(3,msgIndx) = -dstLocalID +! +! end do +! end do +! +! halo%numLocalCopies = msgIndx endif @@ -4633,7 +4646,7 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) case ('northwest') - !*** normal northeast boundary - just copy NW corner + !*** normal northwest boundary - just copy NW corner !*** of physical domain into SE halo of NW nbr block if (dstBlock > 0) then diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 index 231fc9dd3..9355d6164 100644 --- a/cicecore/drivers/unittest/halochk/halochk.F90 +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -5,7 +5,7 @@ module halochk_data use ice_kinds_mod, only: int_kind, dbl_kind, real_kind, log_kind use ice_blocks, only: block, get_block, nx_block, ny_block, nblocks_tot, nghost use ice_boundary, only: ice_HaloUpdate - use ice_constants, only: c0, c1, & + use ice_constants, only: c0, c1, p5, & field_loc_unknown, field_loc_noupdate, & field_loc_center, field_loc_NEcorner, & field_loc_Nface, field_loc_Eface, & @@ -28,6 +28,7 @@ module halochk_data end module halochk_data +!======================================================================= program halochk @@ -39,8 +40,9 @@ program halochk implicit none - integer(int_kind) :: nn, nl, nt, i, j, k1, k2, n, ib, ie, jb, je, iblock - integer(int_kind) :: blockID, numBlocks + integer(int_kind) :: nn, nl, nt, i, j, k1, k2, n, ib, ie, jb, je + integer(int_kind) :: iblock, itrip, ioffset, joffset + integer(int_kind) :: blockID, numBlocks, jtrip type (block) :: this_block real(dbl_kind) , allocatable :: darrayi1(:,:,:) , darrayj1(:,:,:) @@ -63,15 +65,16 @@ program halochk integer(int_kind), parameter :: maxlocs = 5 integer(int_kind), parameter :: nz1 = 3 integer(int_kind), parameter :: nz2 = 4 - real(dbl_kind) :: aichk,ajchk,cichk,cjchk + real(dbl_kind) :: aichk,ajchk,cichk,cjchk,rival,rjval,rsign character(len=16) :: locs_name(maxlocs), types_name(maxtypes) integer(int_kind) :: field_loc(maxlocs), field_type(maxtypes) - integer(int_kind) :: npes, ierr, ntask, testcnt, tottest - integer(int_kind) :: errorflag0, gflag, k1m, k2m, ptcntsum + integer(int_kind) :: npes, ierr, ntask, testcnt, tottest, tpcnt, tfcnt + integer(int_kind) :: errorflag0, gflag, k1m, k2m, ptcntsum, failcntsum integer(int_kind), allocatable :: errorflag(:) - integer(int_kind), allocatable :: ptcnt(:) - character(len=128), allocatable :: stringflag(:) - character(len=32) :: string + integer(int_kind), allocatable :: ptcnt(:), failcnt(:) + character(len=128), allocatable :: teststring(:) + character(len=32) :: halofld + logical :: tripole_average, tripole_pole, spvalL1 logical :: first_call = .true. real(dbl_kind) , parameter :: fillval = -88888.0_dbl_kind @@ -118,9 +121,11 @@ program halochk tottest = maxtests * maxlocs * maxtypes allocate(errorflag(tottest)) - allocate(stringflag(tottest)) + allocate(teststring(tottest)) allocate(ptcnt(tottest)) + allocate(failcnt(tottest)) ptcnt(:) = 0 + failcnt(:) = 0 !----------------------------------------------------------------- ! Testing @@ -145,7 +150,7 @@ program halochk errorflag0 = passflag errorflag(:) = passflag - stringflag(:) = ' ' + teststring(:) = ' ' ! --------------------------- ! TEST HALO UPDATE @@ -295,9 +300,11 @@ program halochk enddo endif - !--- halo off on north boundary --- - if (ns_boundary_type == 'closed' .or. & - ns_boundary_type == 'open' ) then + !--- halo off on north boundary, tripole handled later --- + if (ns_boundary_type == 'closed' .or. & + ns_boundary_type == 'open' .or. & + ns_boundary_type == 'tripole' .or. & + ns_boundary_type == 'tripoleT' ) then do iblock = 1,numBlocks call ice_distributionGetBlockID(distrb_info, iblock, blockID) this_block = get_block(blockID, blockID) @@ -327,8 +334,8 @@ program halochk write(6,*) ' ' write(6,*) 'HALOCHK FAILED' write(6,*) ' ' - call abort_ice(subname//' testcnt > tottest',file=__FILE__,line=__LINE__) endif + call abort_ice(subname//' testcnt > tottest',file=__FILE__,line=__LINE__) endif !--- fill arrays --- @@ -370,25 +377,25 @@ program halochk if (nn == 1) then k1m = 1 k2m = 1 - string = '2DR8' + halofld = '2DR8' call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) elseif (nn == 2) then k1m = nz1 k2m = 1 - string = '3DR8' + halofld = '3DR8' call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) elseif (nn == 3) then k1m = nz1 k2m = nz2 - string = '4DR8' + halofld = '4DR8' call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) elseif (nn == 4) then k1m = 1 k2m = 1 - string = '2DR4' + halofld = '2DR4' rarrayi1 = real(darrayi1,kind=real_kind) rarrayj1 = real(darrayj1,kind=real_kind) call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) @@ -398,7 +405,7 @@ program halochk elseif (nn == 5) then k1m = nz1 k2m = 1 - string = '3DR4' + halofld = '3DR4' rarrayi2 = real(darrayi2,kind=real_kind) rarrayj2 = real(darrayj2,kind=real_kind) call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) @@ -408,7 +415,7 @@ program halochk elseif (nn == 6) then k1m = nz1 k2m = nz2 - string = '4DR4' + halofld = '4DR4' rarrayi3 = real(darrayi3,kind=real_kind) rarrayj3 = real(darrayj3,kind=real_kind) call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) @@ -418,7 +425,7 @@ program halochk elseif (nn == 7) then k1m = 1 k2m = 1 - string = '2DI4' + halofld = '2DI4' iarrayi1 = nint(darrayi1) iarrayj1 = nint(darrayj1) call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) @@ -428,7 +435,7 @@ program halochk elseif (nn == 8) then k1m = nz1 k2m = 1 - string = '3DI4' + halofld = '3DI4' iarrayi2 = nint(darrayi2) iarrayj2 = nint(darrayj2) call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) @@ -438,7 +445,7 @@ program halochk elseif (nn == 9) then k1m = nz1 k2m = nz2 - string = '4DI4' + halofld = '4DI4' iarrayi3 = nint(darrayi3) iarrayj3 = nint(darrayj3) call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) @@ -448,7 +455,7 @@ program halochk elseif (nn == 10) then k1m = 1 k2m = 1 - string = '2DL1' + halofld = '2DL1' larrayi1 = .true. where (darrayi1 == fillval) larrayi1 = .false. larrayj1 = .false. @@ -461,7 +468,7 @@ program halochk where (larrayj1) darrayj1 = c1 endif - write(stringflag(testcnt),'(5a10)') trim(string),trim(locs_name(nl)),trim(types_name(nt)), & + write(teststring(testcnt),'(5a10)') trim(halofld),trim(locs_name(nl)),trim(types_name(nt)), & trim(ew_boundary_type),trim(ns_boundary_type) do iblock = 1,numBlocks @@ -478,15 +485,25 @@ program halochk do i = ib-nghost, ie+nghost do k1 = 1,k1m do k2 = 1,k2m - if (nn == 1 .or. nn == 4 .or. nn == 7 .or. nn == 10) then + tripole_average = .false. + tripole_pole = .false. + spvalL1 = .false. + if (index(halofld,'2D') > 0) then aichk = darrayi1(i,j,iblock) ajchk = darrayj1(i,j,iblock) - elseif (nn == 2 .or. nn == 5 .or. nn == 8) then + elseif (index(halofld,'3D') > 0) then aichk = darrayi2(i,j,k1,iblock) ajchk = darrayj2(i,j,k1,iblock) - elseif (nn == 3 .or. nn == 6 .or. nn == 9) then + elseif (index(halofld,'4D') > 0) then aichk = darrayi3(i,j,k1,k2,iblock) ajchk = darrayj3(i,j,k1,k2,iblock) + else + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK FAILED' + write(6,*) ' ' + endif + call abort_ice(subname//' halofld not matched '//trim(halofld),file=__FILE__,line=__LINE__) endif if (field_loc (nl) == field_loc_noupdate .or. & @@ -496,14 +513,141 @@ program halochk else cichk = cidata_std(i,j,k1,k2,iblock) cjchk = cjdata_std(i,j,k1,k2,iblock) + + !--- tripole on north boundary, need to hardcode --- + !--- tripole and tripoleT slightly different --- + !--- establish special set of points here --- + if ((this_block%j_glob(je) == ny_global) .and. & + ((ns_boundary_type == 'tripole' .and. & + (j > je .or. & + (j == je .and. (field_loc(nl) == field_loc_Nface .or. field_loc(nl) == field_loc_NEcorner)))) .or. & + (ns_boundary_type == 'tripoleT' .and. & + (j >= je)))) then + + ! flip sign for vector/angle + if (field_type(nt) == field_type_vector .or. field_type(nt) == field_type_angle ) then + rsign = -c1 + else + rsign = c1 + endif + + ! for tripole + if (ns_boundary_type == 'tripole') then + + ! compute itrip and jtrip, these are the location where the halo values are defined for i,j + ! for j=je averaging, itrip and jtrip are the 2nd gridpoint associated with averaging + + ! standard center tripole u-fold + itrip = nx_global-this_block%i_glob(i)+1 + jtrip = max(je - (j-je) + 1 , je) + ioffset = 0 + joffset = 0 + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Nface) then + ! need j offset + joffset = -1 + if (j == je) then + tripole_average = .true. + endif + endif + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Eface) then + ! fold plus cell offset + ioffset = -1 + ! CICE treats j=ny_global tripole edge points incorrectly + ! should do edge wraparound and average + ! CICE does not update those points, assumes it's "land" + if (j == je) then + if (this_block%i_glob(i) == nx_global/2) then + tripole_pole = .true. + elseif (this_block%i_glob(i) == nx_global ) then + tripole_pole = .true. + endif + endif + endif + + ! for tripoleT + elseif (ns_boundary_type == 'tripoleT') then + + ! compute itrip and jtrip, these are the location where the halo values are defined for i,j + ! for j=je averaging, itrip and jtrip are the 2nd gridpoint associated with averaging + + ! standard center tripoleT t-fold + itrip = nx_global-this_block%i_glob(i)+2 + jtrip = je - (j-je) + ioffset = 0 + joffset = 0 + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Eface) then + ! fold plus cell offset + ioffset = -1 + endif + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Nface) then + ! need j offset + joffset = -1 + endif + + if (field_loc(nl) == field_loc_Center .or. field_loc(nl) == field_loc_Eface) then + if (j == je) then + tripole_average = .true. + endif + endif + + ! center point poles need to be treated special + if (field_loc(nl) == field_loc_Center) then + if (j == je .and. & + (this_block%i_glob(i) == 1 .or. this_block%i_glob(i) == nx_global/2+1)) then + tripole_pole = .true. + endif + endif + + endif + + itrip = mod(itrip + ioffset + nx_global-1,nx_global)+1 + jtrip = jtrip + joffset + + rival = (real(itrip,kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind) + rjval = (real(this_block%j_glob(jtrip),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind) + + if (index(halofld,'L1') > 0 .and. j == je) then + ! force cichk and cjchk to match on tripole average index, calc not well defined + spvalL1 = .true. + cichk = aichk + cjchk = ajchk + elseif (tripole_pole) then + ! ends of tripole seam not averaged due in CICE to assumption of land + cichk = rsign * cidata_std(i,j,k1,k2,iblock) + cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) + elseif (tripole_average) then + ! tripole average + cichk = p5 * (cidata_std(i,j,k1,k2,iblock) + rsign * rival) + cjchk = p5 * (cjdata_std(i,j,k1,k2,iblock) + rsign * rjval) + else + ! standard tripole fold + cichk = rsign * rival + cjchk = rsign * rjval + endif + +! if (testcnt == 6 .and. j == 61 .and. i < 3) then +! if (testcnt == 186 .and. j == 61 .and. i<4) then +! if (testcnt == 13 .and. j > 61 .and. (i < 3 .or. i > 89)) then +! if (testcnt == 5 .and. j >= 61 .and. (i < 3 .or. i > 90)) then +! write(100+my_task,'(a,5i6,2l3,f6.2,i6)') 'tcx1 ',i,j,iblock,itrip,jtrip,tripole_average,tripole_pole,rsign,this_block%i_glob(i) +! write(100+my_task,'(a,4f12.2)') 'tcx2 ',cidata_std(i,j,k1,k2,iblock),rival,cichk,aichk +! write(100+my_task,'(a,4f12.2)') 'tcx3 ',cjdata_std(i,j,k1,k2,iblock),rjval,cjchk,ajchk +! endif + endif ! tripole or tripoleT endif - if (nn >= 7 .and. nn <= 9) then + if (index(halofld,'I4') > 0) then cichk = real(nint(cichk),kind=dbl_kind) cjchk = real(nint(cjchk),kind=dbl_kind) endif - if (nn == 10) then + if (index(halofld,'L1') > 0 .and. .not.spvalL1) then if (cichk == dhalofillval .or. cichk == fillval) then cichk = c0 else @@ -517,10 +661,10 @@ program halochk endif ptcnt(testcnt) = ptcnt(testcnt) + 1 - call chkresults(aichk,cichk,errorflag(testcnt),testcnt, & - i,j,k1,k2,iblock,first_call,stringflag(testcnt),trim(string)//'_I') - call chkresults(ajchk,cjchk,errorflag(testcnt),testcnt, & - i,j,k1,k2,iblock,first_call,stringflag(testcnt),trim(string)//'_J') + call chkresults(aichk,cichk,errorflag(testcnt),testcnt,failcnt(testcnt), & + i,j,k1,k2,iblock,first_call,teststring(testcnt),trim(halofld)//'_I') + call chkresults(ajchk,cjchk,errorflag(testcnt),testcnt,failcnt(testcnt), & + i,j,k1,k2,iblock,first_call,teststring(testcnt),trim(halofld)//'_J') enddo ! k2 enddo ! k1 enddo ! i @@ -540,6 +684,8 @@ program halochk errorflag(n) = gflag ptcntsum = global_sum(ptcnt(n),distrb_info) ptcnt(n) = ptcntsum + failcntsum = global_sum(failcnt(n),distrb_info) + failcnt(n) = failcntsum enddo errorflag0 = maxval(errorflag(:)) @@ -547,14 +693,21 @@ program halochk write(6,*) ' ' write(6,*) 'HALOCHK COMPLETED SUCCESSFULLY' write(6,*) ' ' + tpcnt = 0 + tfcnt = 0 do n = 1,tottest if (errorflag(n) == passflag) then - write(6,*) 'PASS ',trim(stringflag(n)),ptcnt(n) + tpcnt = tpcnt + 1 + write(6,*) 'PASS ',trim(teststring(n)),ptcnt(n),failcnt(n) else - write(6,*) 'FAIL ',trim(stringflag(n)),ptcnt(n) + tfcnt = tfcnt + 1 + write(6,*) 'FAIL ',trim(teststring(n)),ptcnt(n),failcnt(n) endif enddo write(6,*) ' ' + write(6,*) ' total pass = ',tpcnt + write(6,*) ' total fail = ',tfcnt + write(6,*) ' ' if (errorflag0 == passflag) then write(6,*) 'HALOCHK TEST COMPLETED SUCCESSFULLY' else @@ -576,30 +729,32 @@ end program halochk !======================================================================= - subroutine chkresults(a1,r1,errorflag,testcnt,i,j,k1,k2,iblock,first_call,stringflag,string) + subroutine chkresults(a1,r1,errorflag,testcnt,failcnt,i,j,k1,k2,iblock,first_call,teststring,halofld) use halochk_data implicit none real(dbl_kind) , intent(in) :: a1,r1 - integer(int_kind), intent(inout) :: errorflag + integer(int_kind), intent(inout) :: errorflag, failcnt integer(int_kind), intent(in) :: i,j,k1,k2,iblock,testcnt logical , intent(inout) :: first_call - character(len=*) , intent(in) :: stringflag,string + character(len=*) , intent(in) :: teststring,halofld logical,parameter :: print_always = .false. character(len=*) , parameter :: subname='(chkresults)' if (a1 /= r1 .or. print_always) then errorflag = failflag + failcnt = failcnt + 1 if (first_call) then write(100+my_task,*) ' ' - write(100+my_task,'(a,i4,2a)') '------- TEST = ',testcnt,' ',trim(stringflag) - write(100+my_task,'(a)') ' test task i j k1 k2 iblock expected computed diff' + write(100+my_task,'(a,i4,2a)') '------- TEST = ',testcnt,' ',trim(teststring) + write(100+my_task,*) ' ' + write(100+my_task,'(a)') ' test task i j k1 k2 iblock expected halocomp diff' first_call = .false. endif - write(100+my_task,1001) trim(string),testcnt,my_task,i,j,k1,k2,iblock,r1,a1,r1-a1 + write(100+my_task,1001) trim(halofld),testcnt,my_task,i,j,k1,k2,iblock,r1,a1,r1-a1 endif 1001 format(a8,7i6,3f12.3) diff --git a/configuration/scripts/options/set_nml.bcyclic b/configuration/scripts/options/set_nml.cyclic similarity index 100% rename from configuration/scripts/options/set_nml.bcyclic rename to configuration/scripts/options/set_nml.cyclic diff --git a/configuration/scripts/options/set_nml.bopen b/configuration/scripts/options/set_nml.open similarity index 100% rename from configuration/scripts/options/set_nml.bopen rename to configuration/scripts/options/set_nml.open diff --git a/configuration/scripts/options/set_nml.tripole b/configuration/scripts/options/set_nml.tripole new file mode 100644 index 000000000..7904b8134 --- /dev/null +++ b/configuration/scripts/options/set_nml.tripole @@ -0,0 +1,3 @@ +grid_type = 'tripole' +ew_boundary_type = 'cyclic' +ns_boundary_type = 'tripole' diff --git a/configuration/scripts/options/set_nml.tripolet b/configuration/scripts/options/set_nml.tripolet new file mode 100644 index 000000000..4bb63dc17 --- /dev/null +++ b/configuration/scripts/options/set_nml.tripolet @@ -0,0 +1,3 @@ +grid_type = 'tripole' +ew_boundary_type = 'cyclic' +ns_boundary_type = 'tripoleT' diff --git a/configuration/scripts/tests/unittest_suite.ts b/configuration/scripts/tests/unittest_suite.ts index 319c91aa6..41bbb7d2d 100644 --- a/configuration/scripts/tests/unittest_suite.ts +++ b/configuration/scripts/tests/unittest_suite.ts @@ -1,14 +1,29 @@ # Test Grid PEs Sets BFB-compare -unittest gx3 1x1 helloworld -unittest gx3 1x1 optargs -unittest gx3 1x1 calchk,short -unittest gx3 4x1x25x29x4 sumchk -unittest gx3 1x1x25x29x16 sumchk -unittest tx1 8x1 sumchk -unittest gx3 4x1 bcstchk -unittest gx3 1x1 bcstchk -unittest gx3 8x2 gridavgchk,dwblockall -unittest gx3 12x1 gridavgchk -unittest gx1 28x1 gridavgchk,dwblockall -unittest gx1 16x2 gridavgchk -unittest gbox128 8x2 gridavgchk +unittest gx3 1x1 helloworld +unittest gx3 1x1 optargs +unittest gx3 1x1 calchk,short +unittest gx3 4x1x25x29x4 sumchk +unittest gx3 1x1x25x29x16 sumchk +unittest tx1 8x1 sumchk +unittest gx3 4x1 bcstchk +unittest gx3 1x1 bcstchk +unittest gx3 8x2 g ridavgchk,dwblockall +unittest gx3 12x1 gridavgchk +unittest gx1 28x1 gridavgchk,dwblockall +unittest gx1 16x2 gridavgchk +unittest gbox128 8x2 gridavgchk +unittest gbox80 1x1x10x10x80 halochk,cyclic,debug +unittest gbox80 1x1x24x23x16 halochk +unittest gbox80 1x1x23x24x16 halochk,cyclic +unittest gbox80 1x1x23x23x16 halochk,open +unittest tx1 1x1x90x60x16 halochk,dwblockall +unittest tx1 1x1x90x60x16 halochk,dwblockall,tripolet +unittest tx1 1x1x95x65x16 halochk,dwblockall +unittest tx1 1x1x95x65x16 halochk,dwblockall,tripolet +unittest gx3 4x2 halochk,dwblockall,debug +unittest gx3 8x2x16x12x10 halochk,cyclic,dwblockall +unittest gx3 17x1x16x12x10 halochk,open,dwblockall +unittest tx1 4x2 halochk,dwblockall +unittest tx1 4x2 halochk,dwblockall,tripolet +unittest tx1 4x2x65x45x10 halochk,dwblockall +unittest tx1 4x2x57x43x12 halochk,dwblockall,tripolet From a934b9313ca0c8ab97e8c932e8bbf2cfff5bcf81 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 14:30:29 -0700 Subject: [PATCH 03/13] - Add haloUpdate_stress test to halo unit test - Add tripoleT support to haloUpdate_stress - Add abort check that nx_global is even for tripole grids - Update documentation --- .../infrastructure/comm/mpi/ice_boundary.F90 | 66 ++++++++++++++----- .../comm/serial/ice_boundary.F90 | 66 ++++++++++++++----- cicecore/cicedyn/infrastructure/ice_grid.F90 | 5 ++ cicecore/drivers/unittest/halochk/halochk.F90 | 47 +++++++++++-- configuration/scripts/tests/unittest_suite.ts | 2 +- doc/source/user_guide/ug_implementation.rst | 45 ++++++++++++- 6 files changed, 189 insertions(+), 42 deletions(-) diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 index 58396a58e..2a7d68c11 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 @@ -5474,6 +5474,7 @@ end subroutine ice_HaloUpdate4DI4 !*********************************************************************** ! This routine updates ghost cells for an input array using ! a second array as needed by the stress fields. +! This is just like 2DR8 except no averaging and only on tripole subroutine ice_HaloUpdate_stress(array1, array2, halo, & fieldLoc, fieldKind, & @@ -5711,30 +5712,61 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & call abort_ice(subname//'ERROR: Unknown field kind') end select - select case (fieldLoc) - case (field_loc_center) ! cell center location + if (halo%tripoleTFlag) then + + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = -1 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location - ioffset = 0 - joffset = 0 + ioffset = 0 + joffset = 1 - case (field_loc_NEcorner) ! cell corner location + case (field_loc_Eface) ! cell center location - ioffset = 1 - joffset = 1 + ioffset = 0 + joffset = 0 - case (field_loc_Eface) + case (field_loc_Nface) ! cell corner (velocity) location - ioffset = 1 - joffset = 0 + ioffset = -1 + joffset = 1 - case (field_loc_Nface) + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select - ioffset = 0 - joffset = 1 + else ! tripole u-fold - case default - call abort_ice(subname//'ERROR: Unknown field location') - end select + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = 0 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location + + ioffset = 1 + joffset = 1 + + case (field_loc_Eface) + + ioffset = 1 + joffset = 0 + + case (field_loc_Nface) + + ioffset = 0 + joffset = 1 + + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select + + endif !*** copy out of global tripole buffer into local !*** ghost cells @@ -5765,7 +5797,7 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** out of range and skipped !*** otherwise do the copy - if (jSrc <= nghost+1 .AND. jDst /= -1 ) then + if (jSrc <= halo%tripoleRows .and. jSrc>0 .and. jDst>0) then array1(iDst,jDst,dstBlock) = isign*bufTripoleR8(iSrc,jSrc) endif diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index 30a513430..aaebcfaad 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -3749,6 +3749,7 @@ end subroutine ice_HaloUpdate4DI4 !*********************************************************************** ! This routine updates ghost cells for an input array using ! a second array as needed by the stress fields. +! This is just like 2DR8 except no averaging and only on tripole subroutine ice_HaloUpdate_stress(array1, array2, halo, & fieldLoc, fieldKind, & @@ -3896,30 +3897,61 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & call abort_ice(subname//'ERROR: Unknown field kind') end select - select case (fieldLoc) - case (field_loc_center) ! cell center location + if (halo%tripoleTFlag) then + + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = -1 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location - ioffset = 0 - joffset = 0 + ioffset = 0 + joffset = 1 - case (field_loc_NEcorner) ! cell corner location + case (field_loc_Eface) ! cell center location - ioffset = 1 - joffset = 1 + ioffset = 0 + joffset = 0 - case (field_loc_Eface) + case (field_loc_Nface) ! cell corner (velocity) location - ioffset = 1 - joffset = 0 + ioffset = -1 + joffset = 1 - case (field_loc_Nface) + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select - ioffset = 0 - joffset = 1 + else ! tripole u-fold - case default - call abort_ice(subname//'ERROR: Unknown field location') - end select + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = 0 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location + + ioffset = 1 + joffset = 1 + + case (field_loc_Eface) + + ioffset = 1 + joffset = 0 + + case (field_loc_Nface) + + ioffset = 0 + joffset = 1 + + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select + + endif !*** copy out of global tripole buffer into local !*** ghost cells @@ -3950,7 +3982,7 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** out of range and skipped !*** otherwise do the copy - if (jSrc <= nghost+1) then + if (jSrc <= halo%tripoleRows .and. jSrc>0 .and. jDst>0) then array1(iDst,jDst,dstBlock) = isign*bufTripoleR8(iSrc,jSrc) endif diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 0d56d3400..770ee9ed9 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -360,6 +360,11 @@ subroutine init_grid1 file=__FILE__, line=__LINE__) endif + if (grid_type == 'tripole' .and. (mod(nx_global,2)/=0)) then + call abort_ice(subname//'ERROR: grid_type tripole requires even nx_global number', & + file=__FILE__, line=__LINE__) + endif + if (trim(grid_type) == 'displaced_pole' .or. & trim(grid_type) == 'tripole' .or. & trim(grid_type) == 'regional' ) then diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 index 9355d6164..1d683009e 100644 --- a/cicecore/drivers/unittest/halochk/halochk.F90 +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -4,7 +4,7 @@ module halochk_data use CICE_InitMod use ice_kinds_mod, only: int_kind, dbl_kind, real_kind, log_kind use ice_blocks, only: block, get_block, nx_block, ny_block, nblocks_tot, nghost - use ice_boundary, only: ice_HaloUpdate + use ice_boundary, only: ice_HaloUpdate, ice_HaloUpdate_stress use ice_constants, only: c0, c1, p5, & field_loc_unknown, field_loc_noupdate, & field_loc_center, field_loc_NEcorner, & @@ -55,12 +55,13 @@ program halochk integer(int_kind), allocatable :: iarrayi2(:,:,:,:) , iarrayj2(:,:,:,:) integer(int_kind), allocatable :: iarrayi3(:,:,:,:,:), iarrayj3(:,:,:,:,:) logical(log_kind), allocatable :: larrayi1(:,:,:) , larrayj1(:,:,:) + real(dbl_kind) , allocatable :: darrayi1str(:,:,:) , darrayj1str(:,:,:) real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) real(dbl_kind), allocatable :: cidata_nup(:,:,:,:,:),cjdata_nup(:,:,:,:,:) real(dbl_kind), allocatable :: cidata_std(:,:,:,:,:),cjdata_std(:,:,:,:,:) - integer(int_kind), parameter :: maxtests = 10 + integer(int_kind), parameter :: maxtests = 11 integer(int_kind), parameter :: maxtypes = 4 integer(int_kind), parameter :: maxlocs = 5 integer(int_kind), parameter :: nz1 = 3 @@ -178,6 +179,9 @@ program halochk allocate(iarrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) allocate(larrayi1 (nx_block,ny_block,max_blocks)) allocate(larrayj1 (nx_block,ny_block,max_blocks)) + allocate(darrayi1str(nx_block,ny_block,max_blocks)) + allocate(darrayj1str(nx_block,ny_block,max_blocks)) + allocate(cidata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) allocate(cjdata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) allocate(cidata_std(nx_block,ny_block,nz1,nz2,max_blocks)) @@ -205,6 +209,8 @@ program halochk iarrayj3 = fillval larrayi1 = .false. larrayj1 = .true. + darrayi1str = fillval + darrayj1str = fillval cidata_bas = fillval cjdata_bas = fillval cidata_std = fillval @@ -345,6 +351,8 @@ program halochk darrayj2(:,:,:,:) = fillval darrayi3(:,:,:,:,:) = fillval darrayj3(:,:,:,:,:) = fillval + darrayi1str(:,:,:) = fillval + darrayj1str(:,:,:) = fillval do iblock = 1,numBlocks call ice_distributionGetBlockID(distrb_info, iblock, blockID) this_block = get_block(blockID, blockID) @@ -466,6 +474,14 @@ program halochk where (larrayi1) darrayi1 = c1 darrayj1 = c0 where (larrayj1) darrayj1 = c1 + elseif (nn == 11) then + k1m = 1 + k2m = 1 + halofld = 'STRESS' + darrayi1str = darrayi1 + darrayj1str = darrayj1 + call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) endif write(teststring(testcnt),'(5a10)') trim(halofld),trim(locs_name(nl)),trim(types_name(nt)), & @@ -491,6 +507,9 @@ program halochk if (index(halofld,'2D') > 0) then aichk = darrayi1(i,j,iblock) ajchk = darrayj1(i,j,iblock) + elseif (index(halofld,'STRESS') > 0) then + aichk = darrayi1(i,j,iblock) + ajchk = darrayj1(i,j,iblock) elseif (index(halofld,'3D') > 0) then aichk = darrayi2(i,j,k1,iblock) ajchk = darrayj2(i,j,k1,iblock) @@ -506,6 +525,7 @@ program halochk call abort_ice(subname//' halofld not matched '//trim(halofld),file=__FILE__,line=__LINE__) endif + if (field_loc (nl) == field_loc_noupdate .or. & field_type(nt) == field_type_noupdate) then cichk = cidata_nup(i,j,k1,k2,iblock) @@ -514,6 +534,12 @@ program halochk cichk = cidata_std(i,j,k1,k2,iblock) cjchk = cjdata_std(i,j,k1,k2,iblock) + if (index(halofld,'STRESS') > 0) then + ! only updates on tripole zipper for tripole grids, use daarayi1str as baseline + cichk = darrayi1str(i,j,iblock) + cjchk = darrayj1str(i,j,iblock) + endif + !--- tripole on north boundary, need to hardcode --- !--- tripole and tripoleT slightly different --- !--- establish special set of points here --- @@ -612,13 +638,23 @@ program halochk rjval = (real(this_block%j_glob(jtrip),kind=dbl_kind) + & real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind) - if (index(halofld,'L1') > 0 .and. j == je) then + if (index(halofld,'STRESS') > 0) then + ! only updates on tripole zipper for tripole grids, not tripoleT + if (tripole_pole) then + ! ends of tripole seam not averaged in CICE + cichk = rsign * cidata_std(i,j,k1,k2,iblock) + cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) + else + cichk = rsign * rival + cjchk = rsign * rjval + endif + elseif (index(halofld,'L1') > 0 .and. j == je) then ! force cichk and cjchk to match on tripole average index, calc not well defined spvalL1 = .true. cichk = aichk cjchk = ajchk elseif (tripole_pole) then - ! ends of tripole seam not averaged due in CICE to assumption of land + ! ends of tripole seam not averaged in CICE cichk = rsign * cidata_std(i,j,k1,k2,iblock) cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) elseif (tripole_average) then @@ -635,7 +671,8 @@ program halochk ! if (testcnt == 186 .and. j == 61 .and. i<4) then ! if (testcnt == 13 .and. j > 61 .and. (i < 3 .or. i > 89)) then ! if (testcnt == 5 .and. j >= 61 .and. (i < 3 .or. i > 90)) then -! write(100+my_task,'(a,5i6,2l3,f6.2,i6)') 'tcx1 ',i,j,iblock,itrip,jtrip,tripole_average,tripole_pole,rsign,this_block%i_glob(i) +! write(100+my_task,'(a,5i6,2l3,f6.2,i6)') 'tcx1 ',i,j,iblock,itrip,jtrip, & +! tripole_average,tripole_pole,rsign,this_block%i_glob(i) ! write(100+my_task,'(a,4f12.2)') 'tcx2 ',cidata_std(i,j,k1,k2,iblock),rival,cichk,aichk ! write(100+my_task,'(a,4f12.2)') 'tcx3 ',cjdata_std(i,j,k1,k2,iblock),rjval,cjchk,ajchk ! endif diff --git a/configuration/scripts/tests/unittest_suite.ts b/configuration/scripts/tests/unittest_suite.ts index 41bbb7d2d..7486e87aa 100644 --- a/configuration/scripts/tests/unittest_suite.ts +++ b/configuration/scripts/tests/unittest_suite.ts @@ -7,7 +7,7 @@ unittest gx3 1x1x25x29x16 sumchk unittest tx1 8x1 sumchk unittest gx3 4x1 bcstchk unittest gx3 1x1 bcstchk -unittest gx3 8x2 g ridavgchk,dwblockall +unittest gx3 8x2 gridavgchk,dwblockall unittest gx3 12x1 gridavgchk unittest gx1 28x1 gridavgchk,dwblockall unittest gx1 16x2 gridavgchk diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 5ed2092c0..fe2d081d3 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -284,7 +284,27 @@ fold, on the logical grid. The point with index i along the ghost T-row of index :math:`n+1` physically coincides with point :math:`{\tt nx\_global}-{\tt i}+1` on the T-row of index :math:`n`. The ghost U-row of index :math:`n+1` physically coincides with the U-row of -index :math:`n-1`. +index :math:`n-1`. In the schematics below, symbols A-H represent +grid points from 1:nx_global at a given j index and the setup of the +tripole seam is depicted within a few rows of the seam. + +.. _tab-tripole: + +.. table:: Tripole (u-fold) Grid Schematic + :align: center + + +--------------+---------------------------------------+--------------+ + | global j | | global j | + ! index | grid point IDs (i index) | index source | + +==============+====+====+====+====+====+====+====+====+==============+ + | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | + +==============+====+====+====+====+====+====+====+====+==============+ + | ny_global+1 | H | G | F | E | D | C | B | A | ny_global | + +==============+====+====+====+====+====+====+====+====+==============+ + | ny_global | A | B | C | D | E | F | G | H | | + +==============+====+====+====+====+====+====+====+====+==============+ + | ny_global-1 | A | B | C | D | E | F | G | H | | + +==============+====+====+====+====+====+====+====+====+==============+ In the T-fold tripole grid, the poles have T-index 1 and and :math:`{\tt nx\_global}/2+1` on the top T-row of the physical grid, and @@ -301,6 +321,26 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row :math:`n-1`, and ghost U-row :math:`n+1` coincides with U-row :math:`n-2`. +.. _tab-tripole: + +.. table:: TripoleT (t-fold) Grid Schematic + :align: center + + +--------------+--------------------------------------------+--------------+ + | global j | | global j | + ! index | grid point IDs (i index) | index source | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global+1 | | H | G | F | E | D | C | B | A | ny_global-1 | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global | A | BH | CG | DF | E | FD | GC | HB | | | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global-1 | A | B | C | D | E | F | G | H | | | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global-2 | A | B | C | D | E | F | G | H | | | + +==============+====+====+====+====+====+====+====+====+====+==============+ + The tripole grid thus requires two special kinds of treatment for certain rows, arranged by the halo-update routines. First, within rows along the fold, coincident points must always have the same value. This @@ -310,7 +350,8 @@ the coincident physical rows. Both operations involve the tripole buffer, which is used to assemble the data for the affected rows. Special treatment is also required in the scattering routine, and when computing global sums one of each pair of coincident points has to be -excluded. +excluded. Halos of center, east, north, and northeast points is supported +generally and require unique halo indexing. ***************** Rectangular grids From 86f54d1539c801e8bf7f320bceccfbe14521df82 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 14:32:43 -0700 Subject: [PATCH 04/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index fe2d081d3..0664094a5 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -294,7 +294,7 @@ tripole seam is depicted within a few rows of the seam. :align: center +--------------+---------------------------------------+--------------+ - | global j | | global j | + | global j | | global j | ! index | grid point IDs (i index) | index source | +==============+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | @@ -327,7 +327,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row :align: center +--------------+--------------------------------------------+--------------+ - | global j | | global j | + | global j | | global j | ! index | grid point IDs (i index) | index source | +==============+====+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | From 675adcf02a7396fe0900dca3a2c2ace8e978ac52 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 16:24:27 -0700 Subject: [PATCH 05/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 0664094a5..e59098ffd 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -306,6 +306,7 @@ tripole seam is depicted within a few rows of the seam. | ny_global-1 | A | B | C | D | E | F | G | H | | +==============+====+====+====+====+====+====+====+====+==============+ + In the T-fold tripole grid, the poles have T-index 1 and and :math:`{\tt nx\_global}/2+1` on the top T-row of the physical grid, and points with T-index i and :math:`{\tt nx\_global}-{\tt i}+2` are @@ -321,7 +322,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row :math:`n-1`, and ghost U-row :math:`n+1` coincides with U-row :math:`n-2`. -.. _tab-tripole: +.. _tab-tripoleT: .. table:: TripoleT (t-fold) Grid Schematic :align: center @@ -341,6 +342,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row | ny_global-2 | A | B | C | D | E | F | G | H | | | +==============+====+====+====+====+====+====+====+====+====+==============+ + The tripole grid thus requires two special kinds of treatment for certain rows, arranged by the halo-update routines. First, within rows along the fold, coincident points must always have the same value. This From b027ef2e574f0615a6ffcd326087916d2bf1ab79 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 16:28:04 -0700 Subject: [PATCH 06/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index e59098ffd..39f2d5474 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -298,13 +298,13 @@ tripole seam is depicted within a few rows of the seam. ! index | grid point IDs (i index) | index source | +==============+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | - +==============+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | H | G | F | E | D | C | B | A | ny_global | - +==============+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global | A | B | C | D | E | F | G | H | | - +==============+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global-1 | A | B | C | D | E | F | G | H | | - +==============+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+--------------+ In the T-fold tripole grid, the poles have T-index 1 and and @@ -332,15 +332,15 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row ! index | grid point IDs (i index) | index source | +==============+====+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | | H | G | F | E | D | C | B | A | ny_global-1 | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global | A | BH | CG | DF | E | FD | GC | HB | | | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global-1 | A | B | C | D | E | F | G | H | | | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global-2 | A | B | C | D | E | F | G | H | | | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ The tripole grid thus requires two special kinds of treatment for From 09668e3a6a42576b299483a25bdebd948a8fe5a3 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 16:33:04 -0700 Subject: [PATCH 07/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 39f2d5474..47f93df18 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -296,7 +296,7 @@ tripole seam is depicted within a few rows of the seam. +--------------+---------------------------------------+--------------+ | global j | | global j | ! index | grid point IDs (i index) | index source | - +==============+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | H | G | F | E | D | C | B | A | ny_global | @@ -330,7 +330,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row +--------------+--------------------------------------------+--------------+ | global j | | global j | ! index | grid point IDs (i index) | index source | - +==============+====+====+====+====+====+====+====+====+====+==============+ + +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | | H | G | F | E | D | C | B | A | ny_global-1 | From 26f07b57dccb82d48b5dccdf51360a8084cd7629 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 17:22:41 -0700 Subject: [PATCH 08/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 47f93df18..e6d53a74c 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -295,7 +295,7 @@ tripole seam is depicted within a few rows of the seam. +--------------+---------------------------------------+--------------+ | global j | | global j | - ! index | grid point IDs (i index) | index source | + | index | grid point IDs (i index) | index source | +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | +--------------+----+----+----+----+----+----+----+----+--------------+ @@ -329,7 +329,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row +--------------+--------------------------------------------+--------------+ | global j | | global j | - ! index | grid point IDs (i index) | index source | + | index | grid point IDs (i index) | index source | +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | +--------------+----+----+----+----+----+----+----+----+----+--------------+ From aa8218cd3288a4beabec58f01f60de746108dbf6 Mon Sep 17 00:00:00 2001 From: apcraig Date: Fri, 10 Mar 2023 17:33:24 -0700 Subject: [PATCH 09/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index e6d53a74c..1a89b90db 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -296,7 +296,7 @@ tripole seam is depicted within a few rows of the seam. +--------------+---------------------------------------+--------------+ | global j | | global j | | index | grid point IDs (i index) | index source | - +--------------+----+----+----+----+----+----+----+----+--------------+ + +==============+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | +--------------+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | H | G | F | E | D | C | B | A | ny_global | @@ -330,7 +330,7 @@ north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row +--------------+--------------------------------------------+--------------+ | global j | | global j | | index | grid point IDs (i index) | index source | - +--------------+----+----+----+----+----+----+----+----+----+--------------+ + +==============+====+====+====+====+====+====+====+====+====+==============+ | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | +--------------+----+----+----+----+----+----+----+----+----+--------------+ | ny_global+1 | | H | G | F | E | D | C | B | A | ny_global-1 | From 611e4001f9df437fb6174b06f0cb70c7650b247b Mon Sep 17 00:00:00 2001 From: apcraig Date: Sat, 11 Mar 2023 13:20:39 -0700 Subject: [PATCH 10/13] Update halochk test to make haloupdate_stress test more robust, less chance for false positive --- cicecore/drivers/unittest/halochk/halochk.F90 | 29 +++++++++++++------ 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 index 1d683009e..29eaa8150 100644 --- a/cicecore/drivers/unittest/halochk/halochk.F90 +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -56,6 +56,7 @@ program halochk integer(int_kind), allocatable :: iarrayi3(:,:,:,:,:), iarrayj3(:,:,:,:,:) logical(log_kind), allocatable :: larrayi1(:,:,:) , larrayj1(:,:,:) real(dbl_kind) , allocatable :: darrayi1str(:,:,:) , darrayj1str(:,:,:) + real(dbl_kind) , allocatable :: darrayi10(:,:,:) , darrayj10(:,:,:) real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) real(dbl_kind), allocatable :: cidata_nup(:,:,:,:,:),cjdata_nup(:,:,:,:,:) @@ -181,6 +182,8 @@ program halochk allocate(larrayj1 (nx_block,ny_block,max_blocks)) allocate(darrayi1str(nx_block,ny_block,max_blocks)) allocate(darrayj1str(nx_block,ny_block,max_blocks)) + allocate(darrayi10 (nx_block,ny_block,max_blocks)) + allocate(darrayj10 (nx_block,ny_block,max_blocks)) allocate(cidata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) allocate(cjdata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) @@ -211,6 +214,8 @@ program halochk larrayj1 = .true. darrayi1str = fillval darrayj1str = fillval + darrayi10 = fillval + darrayj10 = fillval cidata_bas = fillval cjdata_bas = fillval cidata_std = fillval @@ -379,6 +384,10 @@ program halochk enddo enddo + ! copy original darray1 for "stress" compare + darrayi10 = darrayi1 + darrayj10 = darrayj1 + !--- halo update --- @@ -478,8 +487,8 @@ program halochk k1m = 1 k2m = 1 halofld = 'STRESS' - darrayi1str = darrayi1 - darrayj1str = darrayj1 + darrayi1str = -darrayi1 ! flip sign for testing + darrayj1str = -darrayj1 call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) endif @@ -535,9 +544,10 @@ program halochk cjchk = cjdata_std(i,j,k1,k2,iblock) if (index(halofld,'STRESS') > 0) then - ! only updates on tripole zipper for tripole grids, use daarayi1str as baseline - cichk = darrayi1str(i,j,iblock) - cjchk = darrayj1str(i,j,iblock) + ! only updates on tripole zipper for tripole grids + ! darrayi10 is copy of darrayi1 before halo call + cichk = darrayi10(i,j,iblock) + cjchk = darrayj10(i,j,iblock) endif !--- tripole on north boundary, need to hardcode --- @@ -641,12 +651,13 @@ program halochk if (index(halofld,'STRESS') > 0) then ! only updates on tripole zipper for tripole grids, not tripoleT if (tripole_pole) then + ! flip sign due to sign of darrayi1str ! ends of tripole seam not averaged in CICE - cichk = rsign * cidata_std(i,j,k1,k2,iblock) - cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) + cichk = -rsign * cidata_std(i,j,k1,k2,iblock) + cjchk = -rsign * cjdata_std(i,j,k1,k2,iblock) else - cichk = rsign * rival - cjchk = rsign * rjval + cichk = -rsign * rival + cjchk = -rsign * rjval endif elseif (index(halofld,'L1') > 0 .and. j == je) then ! force cichk and cjchk to match on tripole average index, calc not well defined From 9e8463c71f32f1a2d297b1f176ba9a2760cab3d8 Mon Sep 17 00:00:00 2001 From: apcraig Date: Sat, 11 Mar 2023 13:24:35 -0700 Subject: [PATCH 11/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 1a89b90db..06f968e0b 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -352,8 +352,8 @@ the coincident physical rows. Both operations involve the tripole buffer, which is used to assemble the data for the affected rows. Special treatment is also required in the scattering routine, and when computing global sums one of each pair of coincident points has to be -excluded. Halos of center, east, north, and northeast points is supported -generally and require unique halo indexing. +excluded. Halos of center, east, north, and northeast points are supported, +and each requires slightly different halo indexing across the tripole seam. ***************** Rectangular grids From d07d9ef69ee3388443c452f2b1c26c362f9c8160 Mon Sep 17 00:00:00 2001 From: apcraig Date: Sun, 12 Mar 2023 14:27:01 -0600 Subject: [PATCH 12/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 06f968e0b..805089631 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -275,14 +275,14 @@ the namelist variable ``ns_boundary_type``, ‘tripole’ for the U-fold and ‘tripoleT’ for the T-fold grid. In the U-fold tripole grid, the poles have U-index -:math:`{\tt nx\_global}/2` and ``nx_global`` on the top U-row of the -physical grid, and points with U-index i and :math:`{\tt nx\_global-i}` +:math:`nx\_global/2` and :math:`nx_global` on the top U-row of the +physical grid, and points with U-index :math:`i` and :math:`nx\_global-i` are coincident. Let the fold have U-row index :math:`n` on the global grid; this will also be the T-row index of the T-row to the south of the fold. There are ghost (halo) T- and U-rows to the north, beyond the fold, on the logical grid. The point with index i along the ghost T-row of index :math:`n+1` physically coincides with point -:math:`{\tt nx\_global}-{\tt i}+1` on the T-row of index :math:`n`. The +:math:`nx\_global-i+1` on the T-row of index :math:`n`. The ghost U-row of index :math:`n+1` physically coincides with the U-row of index :math:`n-1`. In the schematics below, symbols A-H represent grid points from 1:nx_global at a given j index and the setup of the @@ -307,16 +307,16 @@ tripole seam is depicted within a few rows of the seam. +--------------+----+----+----+----+----+----+----+----+--------------+ -In the T-fold tripole grid, the poles have T-index 1 and and -:math:`{\tt nx\_global}/2+1` on the top T-row of the physical grid, and -points with T-index i and :math:`{\tt nx\_global}-{\tt i}+2` are +In the T-fold tripole grid, the poles have T-index :math:`1` and and +:math:`nx\_global/2+1` on the top T-row of the physical grid, and +points with T-index :math:`i` and :math:`nx\_global-i+2` are coincident. Let the fold have T-row index :math:`n` on the global grid. It is usual for the northernmost row of the physical domain to be a U-row, but in the case of the T-fold, the U-row of index :math:`n` is “beyond” the fold; although it is not a ghost row, it is not physically independent, because it coincides with U-row :math:`n-1`, and it therefore has to be treated like a ghost row. Points i on U-row -:math:`n` coincides with :math:`{\tt nx\_global}-{\tt i}+1` on U-row +:math:`n` coincides with :math:`nx\_global-i+1` on U-row :math:`n-1`. There are still ghost T- and U-rows :math:`n+1` to the north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row :math:`n-1`, and ghost U-row :math:`n+1` coincides with U-row From 4bba4575dddb9ca730cd72f9f8dec1fa94049718 Mon Sep 17 00:00:00 2001 From: apcraig Date: Sun, 12 Mar 2023 14:30:09 -0600 Subject: [PATCH 13/13] update documentation --- doc/source/user_guide/ug_implementation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 805089631..e79130e02 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -275,7 +275,7 @@ the namelist variable ``ns_boundary_type``, ‘tripole’ for the U-fold and ‘tripoleT’ for the T-fold grid. In the U-fold tripole grid, the poles have U-index -:math:`nx\_global/2` and :math:`nx_global` on the top U-row of the +:math:`nx\_global/2` and :math:`nx\_global` on the top U-row of the physical grid, and points with U-index :math:`i` and :math:`nx\_global-i` are coincident. Let the fold have U-row index :math:`n` on the global grid; this will also be the T-row index of the T-row to the south of the