From a2dafd25d2c9c81e08fafddc6f1f20f5c9b0dd5f Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Thu, 5 Aug 2021 09:58:29 -0600 Subject: [PATCH] Changes to support CESM ocean coupling This commit includes several changes to support CISM coupling with POP in CESM. * Modified the ocean_data_domain option, which determines where thermal forcing is computed when running with whichbmlt_float = BMLT_FLOAT_THERMAL_FORCING = 6. There are still three values, but they have new meanings: (0) ocean data computed internally by CISM (1) ocean data read from external file (2) ocean data received from coupler via Glad Option (0) is not yet supported, but would be used with an internal plume model. Option (1) is the default and is used for ISMIP6 Antarctic runs. Option (2) should be used for CESM coupled runs. * Added an ocean_data_extrapolate option: (0) Do not extrapolate ocean thermal forcing to ice-shelf cavities (1) Extrapolate ocean thermal forcing to ice-shelf cavities Option (0) is the default and is used for ISMIP6 runs where thermal forcing is provided for the entire Antarctic domains. Since it is the default, existing ISMIP6 config files do not need to be modified. Option (1) should be used for CESM coupled runs in which the ocean does not supply thermal forcing in ice-shelf cavities, as will always be the case with POP. * Modified logic in subroutine glad_i_tstep_gcm of glad_timestep.F90 so that TF from Glad does not overwrite TF from the CISM input file unless ocean_data_domain = 2. * Added a parallel interface, parallel_boundary_value, that inserts a user-specified value in grid cells at the global boundary. This subroutine is used to set TF = unphys_val at the global boundary after halo calls for TF in module glissade_bmlt_float. The rationale is as follows: CISM Antarctic runs often use global no_ice BCs to enable running on active blocks only. These BCs require that scalars are set to 0 in the row of cells at the global boundary, just inside the global halo. If TF = 0 in these cells, but TF = unphys_val in all other cells, the extrapolation routine will set TF = 0 in neighboring cells until putting apparently realistic (but actually fictitious) values of TF throughout the domain. It would be better *not* to put 0 values in boundary cells with the no_ice BCs; might want to return to this later. * Added a warning message if whichbmlt_float = 6, but TF <= 0 everywhere. * Removed code that aborts CISM when the extrapolated TF < 0. We have TF <= 0 everywhere for ISMIP6 datasets, but this will not be true in general. * Added code that aborts CISM whether TF has an extreme value, outside the range (-5K, 20K). * Fixed a diagnostic print that was trying to write the k index when k was undefined. * Cleaned up some other comments and diagnostics. With these changes, the code runs to completion in cism-wrapper tests with a dynamic Antarctic ice sheet, both when reading TF from an input file and when deriving TF from idealized T and S fields set in Glad. It remains to test the code when receiving T and S from POP. Note: To run cism-wrapper tests, I modified the ocean_data_domain entry and added an ocean_data_extrapolate entry in .../cime_config/namelist_definition_cism.xml. I confirmed that these changes are BFB for a typical ISMIP6 Antarctic run on Cheyenne. --- libglad/glad_main.F90 | 21 +++-- libglad/glad_mbal_coupling.F90 | 8 +- libglad/glad_timestep.F90 | 20 ++-- libglide/glide_setup.F90 | 14 ++- libglide/glide_types.F90 | 23 +++-- libglimmer/parallel_mpi.F90 | 82 +++++++++++++++- libglimmer/parallel_slap.F90 | 86 ++++++++++++++++- libglissade/glissade.F90 | 29 +++++- libglissade/glissade_bmlt_float.F90 | 139 +++++++++++++--------------- 9 files changed, 307 insertions(+), 115 deletions(-) diff --git a/libglad/glad_main.F90 b/libglad/glad_main.F90 index b2e19a44..d2fae8d0 100644 --- a/libglad/glad_main.F90 +++ b/libglad/glad_main.F90 @@ -125,10 +125,10 @@ module glad_main ! When coupled to CESM, Glad receives several fields from the coupler on the ice sheet grid: ! qsmb = surface mass balance (kg/m^2/s) ! tsfc = surface ground temperature (deg C) - ! salinity1..7 = ocean salinity at levels 0,10,19,26,30,33,35 (g/kg) - ! tocn1..7 = ocean temperatures at levels 0,10,19,26,30,33,35 (deg K) + ! salinity = ocean salinity at one or more ocean levels (g/kg) + ! tocn = ocean temperature at one or more levels (deg K) ! Both qsmb and tsfc are computed in the CESM land model. - ! Both set of fields salinity1..7 and tocan1..7 are computed in POP. + ! Both salinity and tocn are computed in the CESM ocean model. ! Seven fields are returned to CESM on the ice sheet grid: ! ice_covered = whether a grid cell is ice-covered [0,1] ! topo = surface elevation (m) @@ -584,13 +584,17 @@ subroutine glad_gcm(params, instance_index, time, & ! Main Glad subroutine for GCM coupling. ! - ! It does all necessary temporal averaging, - ! and calls the dynamic ice sheet model when required. + ! It does all temporal averaging and calls the dynamic ice sheet model when required. ! ! Input fields should be taken as means over the period since the last call. ! See the user documentation for more information. ! - ! Input fields are assumed to NOT have halo cells + ! Input fields are assumed to NOT have halo cells. + ! + ! Glad accumulates and averages ocean thermal_forcing based on salinity and tocn, + ! but this thermal forcing is used by CISM only if model%options%ocean_data_domain = 2. + ! Otherwise, CISM will compute thermal_forcing internally (ocean_data_domain = 0) + ! or read it from an external file (ocean_data_domain = 1). use glimmer_utils use glad_timestep, only : glad_i_tstep_gcm @@ -711,13 +715,14 @@ subroutine glad_gcm(params, instance_index, time, & thermal_forcing_haloed(k,:,:)) enddo + if (verbose_glad .and. this_rank == rtest) then i = itest j = jtest print*, 'r, i, j =', this_rank, i, j print*, 'k, zocn, temperature, salinity, thermal forcing:' do k = 1, nzocn - write(6,'(i4, 4f10.3)') k, zocn(k), & + write(6,'(i4, 4f11.3)') k, zocn(k), & tocn_haloed(k,i,j), salinity_haloed(k,i,j), thermal_forcing_haloed(k,i,j) enddo endif @@ -799,7 +804,7 @@ subroutine glad_gcm(params, instance_index, time, & j = jtest print*, 'Before calling glad_i_tstep_gcm, k, zocn, average thermal forcing:' do k = 1, nzocn - write(6,'(i4, 2f10.3)') k, zocn(k), params%instances(instance_index)%thermal_forcing(k,i,j) + write(6,'(i4, 2f11.3)') k, zocn(k), params%instances(instance_index)%thermal_forcing(k,i,j) enddo endif diff --git a/libglad/glad_mbal_coupling.F90 b/libglad/glad_mbal_coupling.F90 index ed69ed83..8894f0d5 100644 --- a/libglad/glad_mbal_coupling.F90 +++ b/libglad/glad_mbal_coupling.F90 @@ -104,9 +104,7 @@ end subroutine glad_mbc_init !++++++++++++++++++++++++++++++++++++++++++++++++++++++ - subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing) !, & -! thermal_forcing2, thermal_forcing3, thermal_forcing4, & -! thermal_forcing5, thermal_forcing6, thermal_forcing7) + subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing) ! In glint, this was done in glint_downscale.F90 @@ -115,7 +113,7 @@ subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing) real(dp),dimension(:,:),intent(in) :: acab ! Surface mass balance (m) real(dp),dimension(:,:),intent(in) :: artm ! Mean air temperature (degC) - real(dp),dimension(:,:,:),intent(in) :: thermal_forcing ! Mean thermal_forcing at level 0 (degK) + real(dp),dimension(:,:,:),intent(in) :: thermal_forcing ! Ocean thermal_forcing (degK) ! Things to do the first time @@ -161,7 +159,7 @@ subroutine glad_average_input_gcm(params, dt, acab, artm, thermal_forcing) integer, intent(in) :: dt !> mbal accumulation time (hours) real(dp),dimension(:,:),intent(out) :: artm !> Mean air temperature (degC) real(dp),dimension(:,:),intent(out) :: acab !> Mass-balance (m/yr) - real(dp),dimension(:,:,:),intent(out) :: thermal_forcing ! Mean thermal_forcing at level 0 (degK) + real(dp),dimension(:,:,:),intent(out) :: thermal_forcing ! Ocean thermal_forcing (degK) if (.not. params%new_accum) then params%artm_save = params%artm_save / real(params%av_count,dp) diff --git a/libglad/glad_timestep.F90 b/libglad/glad_timestep.F90 index 97e66783..7119d262 100644 --- a/libglad/glad_timestep.F90 +++ b/libglad/glad_timestep.F90 @@ -66,6 +66,7 @@ subroutine glad_i_tstep_gcm(time, instance, & use glide use glissade use glide_io + use glide_types use glad_mbal_coupling, only : glad_accumulate_input_gcm, glad_average_input_gcm use glad_io use glad_mbal_io @@ -200,16 +201,21 @@ subroutine glad_i_tstep_gcm(time, instance, & call glide_set_acab(instance%model, instance%acab * rhow/rhoi) call glide_set_artm(instance%model, instance%artm) + ! Note: CISM has several thermal forcing options, as determined by ocean_data_domain: + ! compute internally, read from external file, or receive from Glad. + ! Only if receiving from Glad do we set instance%model%ocean_data%thermal_forcing here. + ! ! Note: The ocean thermal forcing is reset only on the first ice dynamics timestep within this loop. - ! The reason is that CISM has the option of extrapolating thermal forcing from open ocean - ! (i.e., the overlap region between the ocean and ice sheet domains) into sub-shelf cavities. - ! It is more efficient to do this extrapolation only once within the mass_balance timestep - ! (with minor subsequent adjustments if CISM ice shelves retreat) than to extrapolate - ! from the open ocean every ice dynamics time step. + ! The reason is that CISM has the option of extrapolating thermal forcing from open ocean + ! (i.e., the overlap region between the ocean and ice sheet domains) into sub-shelf cavities. + ! It is more efficient to do this extrapolation only once within the mass_balance timestep + ! (with minor subsequent adjustments if CISM ice shelves retreat) than to extrapolate + ! from the open ocean every ice dynamics time step. if (iter == 1) then - if ( associated(instance%model%ocean_data%thermal_forcing) ) then - ! GL: At this point, glide_set does not work for 3D variables. + if ( associated(instance%model%ocean_data%thermal_forcing) .and. & + instance%model%options%ocean_data_domain == OCEAN_DATA_GLAD) then + ! GRL: At this point, glide_set does not work for 3D variables. instance%model%ocean_data%thermal_forcing = instance%thermal_forcing endif endif diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index e087d2eb..3b967d05 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -722,6 +722,7 @@ subroutine handle_options(section, model) call GetValue(section,'bmlt_float_thermal_forcing_param',model%options%bmlt_float_thermal_forcing_param) call GetValue(section,'bmlt_float_ismip6_magnitude',model%options%bmlt_float_ismip6_magnitude) call GetValue(section,'ocean_data_domain',model%options%ocean_data_domain) + call GetValue(section,'ocean_data_extrapolate',model%options%ocean_data_extrapolate) call GetValue(section,'enable_bmlt_anomaly',model%options%enable_bmlt_anomaly) call GetValue(section,'basal_mass_balance',model%options%basal_mbal) call GetValue(section,'smb_input',model%options%smb_input) @@ -936,9 +937,13 @@ subroutine print_options(model) 'highest forcing magnitude ' /) character(len=*), dimension(0:2), parameter :: ocean_data_domain = (/ & - 'data from external ocean model domain; extrapolate to cavities', & - 'ocean data already extrapolated to ice shelf cavities ', & - 'apply data in CISM ice-free ocean; extrapolate to cavities ' /) + 'ocean data computed internally by CISM', & + 'ocean data read from external file ', & + 'ocean data from coupler via Glad ' /) + + character(len=*), dimension(0:1), parameter :: ocean_data_extrapolate = (/ & + 'ocean data not extrapolated to cavities', & + 'ocean data extrapolated to cavities ' /) character(len=*), dimension(0:1), parameter :: smb_input = (/ & 'SMB input in units of m/yr ice ', & @@ -1522,6 +1527,9 @@ subroutine print_options(model) write(message,*) 'ocean data domain : ', model%options%ocean_data_domain, & ocean_data_domain(model%options%ocean_data_domain) call write_log(message) + write(message,*) 'ocean data extrapolate : ', model%options%ocean_data_extrapolate, & + ocean_data_extrapolate(model%options%ocean_data_extrapolate) + call write_log(message) endif if (model%options%basal_mbal < 0 .or. model%options%basal_mbal >= size(b_mbal)) then diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index c633d0d1..6954b5a9 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -137,9 +137,12 @@ module glide_types integer, parameter :: BMLT_FLOAT_ISMIP6_MEDIAN = 1 integer, parameter :: BMLT_FLOAT_ISMIP6_PCT95 = 2 - integer, parameter :: DATA_OCEAN_ONLY = 0 - integer, parameter :: DATA_OCEAN_ICE = 1 - integer, parameter :: DATA_CISM_OCEAN_MASK = 2 + integer, parameter :: OCEAN_DATA_INTERNAL = 0 + integer, parameter :: OCEAN_DATA_EXTERNAL = 1 + integer, parameter :: OCEAN_DATA_GLAD = 2 + + integer, parameter :: OCEAN_DATA_EXTRAPOLATE_FALSE = 0 + integer, parameter :: OCEAN_DATA_EXTRAPOLATE_TRUE = 1 integer, parameter :: BASAL_MBAL_NO_CONTINUITY = 0 integer, parameter :: BASAL_MBAL_CONTINUITY = 1 @@ -522,12 +525,18 @@ module glide_types !> \item[2] High level of forcing (e.g., pct95) !> \end{description} - integer :: ocean_data_domain = 0 + integer :: ocean_data_domain = 1 + + !> \begin{description} + !> \item[0] ocean data computed internally by CISM + !> \item[1] ocean data read from external file + !> \item[2] ocean data received from coupler via Glad + !> \end{description} + integer :: ocean_data_extrapolate = 0 !> \begin{description} - !> \item[0] ocean data on ocean domain only; extrapolate data to shelf cavities - !> \item[1] ocean data available everywhere; already extrapolated to shelf cavities - !> \item[2] ocean data applied where CISM has ice-free ocean; extrapolated to shelf cavities + !> \item[0] ocean data not extrapolated to shelf cavities + !> \item[1] ocean data extrapolated to shelf cavities !> \end{description} logical :: enable_bmlt_anomaly = .false. diff --git a/libglimmer/parallel_mpi.F90 b/libglimmer/parallel_mpi.F90 index 4da72501..fa8ed875 100644 --- a/libglimmer/parallel_mpi.F90 +++ b/libglimmer/parallel_mpi.F90 @@ -248,6 +248,11 @@ module cism_parallel module procedure distributed_scatter_var_col_real8_2d end interface + interface parallel_boundary_value + module procedure parallel_boundary_value_real8_2d + module procedure parallel_boundary_value_real8_3d + end interface parallel_boundary_value + interface parallel_convert_haloed_to_nonhaloed module procedure parallel_convert_haloed_to_nonhaloed_real4_2d module procedure parallel_convert_haloed_to_nonhaloed_real8_2d @@ -5020,7 +5025,7 @@ end subroutine parallel_barrier !======================================================================= - function parallel_boundary(ew, ewn, ns, nsn, parallel) + function parallel_boundary(ew, ns, parallel) implicit none logical :: parallel_boundary @@ -5029,6 +5034,8 @@ function parallel_boundary(ew, ewn, ns, nsn, parallel) ! begin associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn, & global_ewn => parallel%global_ewn, & global_nsn => parallel%global_nsn, & ewlb => parallel%ewlb, & @@ -5038,14 +5045,83 @@ function parallel_boundary(ew, ewn, ns, nsn, parallel) ) parallel_boundary = (ewlb<1 .and. ew==1+lhalo) .or. & - (ewub>global_ewn .and. ew==ewn-uhalo) .or. & + (ewub>global_ewn .and. ew==local_ewn-uhalo) .or. & (nslb<1 .and. ns==1+lhalo) .or. & - (nsub>global_nsn .and. ns==nsn-uhalo) + (nsub>global_nsn .and. ns==local_nsn-uhalo) end associate end function parallel_boundary +!======================================================================= + + ! subroutines belonging to the parallel_boundary_value interface + + subroutine parallel_boundary_value_real8_2d(& + field, & + boundary_value, & + parallel) + + ! Insert a specified value into cells on the global boundary. + ! Typically called with value = 0.0 or an special value. + + real(dp), dimension(:,:), intent(inout) :: field + real(dp), intent(in) :: boundary_value + type(parallel_type) :: parallel + + integer :: ew, ns + + ! begin + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn & + ) + + do ns = 1, local_nsn + do ew = 1, local_ewn + if (parallel_boundary(ew, ns, parallel)) then + field(ew,ns) = boundary_value + endif + enddo + enddo + + end associate + + end subroutine parallel_boundary_value_real8_2d + + + subroutine parallel_boundary_value_real8_3d(& + field, & + boundary_value, & + parallel) + + ! Insert a specified value into cells on the global boundary. + ! Typically called with value = 0.0 or an special value. + + real(dp), dimension(:,:,:), intent(inout) :: field + real(dp), intent(in) :: boundary_value + type(parallel_type) :: parallel + + integer :: ew, ns + + ! begin + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn & + ) + + do ns = 1, local_nsn + do ew = 1, local_ewn + if (parallel_boundary(ew, ns, parallel)) then + field(:,ew,ns) = boundary_value + endif + enddo + enddo + + end associate + + end subroutine parallel_boundary_value_real8_3d + !======================================================================= function parallel_close(ncid) diff --git a/libglimmer/parallel_slap.F90 b/libglimmer/parallel_slap.F90 index b0b88d59..df5991cf 100644 --- a/libglimmer/parallel_slap.F90 +++ b/libglimmer/parallel_slap.F90 @@ -215,6 +215,11 @@ module cism_parallel module procedure distributed_scatter_var_col_real8_2d end interface + interface parallel_boundary_value + module procedure parallel_boundary_value_real8_2d + module procedure parallel_boundary_value_real8_3d + end interface parallel_boundary_value + interface parallel_convert_haloed_to_nonhaloed module procedure parallel_convert_haloed_to_nonhaloed_real4_2d module procedure parallel_convert_haloed_to_nonhaloed_real8_2d @@ -1836,17 +1841,92 @@ end subroutine parallel_barrier !======================================================================= - function parallel_boundary(ew,ewn,ns,nsn) + function parallel_boundary(ew, ns, parallel) implicit none logical :: parallel_boundary - integer :: ew,ewn,ns,nsn + integer :: ew, ns + type(parallel_type) :: parallel + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) ! begin - parallel_boundary = (ew==1 .or. ew==ewn .or. ns==1 .or. ns==nsn) + parallel_boundary = (ew==1 .or. ew==local_ewn .or. & + ns==1 .or. ns==local_nsn) + + end associate end function parallel_boundary +!======================================================================= + + ! subroutines belonging to the parallel_boundary_value interface + + subroutine parallel_boundary_value_real8_2d(& + field, & + boundary_value, & + parallel) + + ! Insert a specified value into cells on the global boundary. + ! Typically called with value = 0.0 or an special value. + + real(dp), dimension(:,:), intent(inout) :: field + real(dp), intent(in) :: boundary_value + type(parallel_type) :: parallel + + integer :: ew, ns + + ! begin + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + do ns = 1, local_nsn + do ew = 1, local_ewn + if (parallel_boundary(ew, ns, parallel)) then + field(ew,ns) = boundary_value + endif + enddo + enddo + + end associate + + end subroutine parallel_boundary_value_real8_2d + + + subroutine parallel_boundary_value_real8_3d(& + field, & + boundary_value, & + parallel) + + ! Insert a specified value into cells on the global boundary. + ! Typically called with value = 0.0 or an special value. + + real(dp), dimension(:,:,:), intent(inout) :: field + real(dp), intent(in) :: boundary_value + type(parallel_type) :: parallel + + integer :: ew, ns + + ! begin + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + do ns = 1, local_nsn + do ew = 1, local_ewn + if (parallel_boundary(ew, ns, parallel)) then + field(:,ew,ns) = boundary_value + endif + enddo + enddo + + end associate + + end subroutine parallel_boundary_value_real8_3d + !======================================================================= function parallel_close(ncid) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index afe4c246..ef0dc5dd 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -1289,12 +1289,13 @@ subroutine glissade_bmlt_float_solve(model) ! Solve for basal melting beneath floating ice. - use glimmer_paramets, only: eps08, tim0, thk0, len0 + use glimmer_paramets, only: eps08, eps11, tim0, thk0, len0 use glimmer_physcon, only: scyr use glissade_bmlt_float, only: glissade_basal_melting_float, & glissade_bmlt_float_thermal_forcing, verbose_bmlt_float use glissade_transport, only: glissade_add_2d_anomaly use glissade_masks, only: glissade_get_masks + use cism_parallel, only: parallel_reduce_max implicit none @@ -1323,6 +1324,7 @@ subroutine glissade_bmlt_float_solve(model) real(dp) :: tf_anomaly_basin ! basin number where anomaly is applied; ! for default value of 0, apply to all basins + real(dp) :: local_maxval, global_maxval ! max values of a given variable integer :: i, j integer :: ewn, nsn real(dp) :: dew, dns @@ -1391,6 +1393,25 @@ subroutine glissade_bmlt_float_solve(model) print*, 'Compute bmlt_float at runtime from current thermal forcing' endif + !Note: Currently, there is no difference between ocean_data_domain = 0 + ! (compute internally) and ocean_data_domain = 1 (read from file). + ! Thermal forcing is initialized to zero and then is loaded from + ! the input or forcing file, if present. + ! If ocean_data_domain = 2, then the thermal forcing is set by Glad; + ! any values read from an input or forcing file are overwritten. + ! CISM is not yet able to compute thermal forcing internally. + !TODO: Add code to compute thermal forcing internally. + + ! Check for positive values of thermal forcing. + ! If whichbmlt_float = BMLT_FLOAT_THERMAL_FORCING, but there are no positive values, + ! something is probably wrong. + + local_maxval = maxval(model%ocean_data%thermal_forcing) + global_maxval = parallel_reduce_max(local_maxval) + if (global_maxval <= eps11) then + call write_log('thermal forcing <= 0 everywhere, GM_WARNING') + endif + !----------------------------------------------- ! Optionally, apply a uniform thermal forcing anomaly everywhere. ! This anomaly can be phased in linearly over a prescribed timescale. @@ -1421,7 +1442,7 @@ subroutine glissade_bmlt_float_solve(model) call glissade_bmlt_float_thermal_forcing(& model%options%bmlt_float_thermal_forcing_param, & - model%options%ocean_data_domain, & + model%options%ocean_data_extrapolate, & parallel, & ewn, nsn, & dew*len0, dns*len0, & ! m @@ -1449,6 +1470,7 @@ subroutine glissade_bmlt_float_solve(model) ! This means that the anomaly is potentially much larger for new cavities ! than for cavities initially present. ! Note: bmlt_float is a basal melting potential; it is reduced below for partly or fully grounded ice. + ! TODO: Remove option (2), which was used for ISMIP6 Antarctica but is now deprecated. if (model%options%which_ho_bmlt_inversion == HO_BMLT_INVERSION_APPLY) then @@ -1492,7 +1514,8 @@ subroutine glissade_bmlt_float_solve(model) ! Convert bmlt_float from SI units (m/s) to scaled model units model%basal_melt%bmlt_float(:,:) = model%basal_melt%bmlt_float(:,:) * tim0/thk0 - else ! other options include BMLT_FLOAT_CONSTANT, BMLT_FLOAT_MISMIP, BMLT_FLOAT_DEPTH, BMLT_FLOAT_MISOMIP and BMLT_FLOAT_POP_CPL + else ! other options include BMLT_FLOAT_CONSTANT, BMLT_FLOAT_MISMIP, & + ! BMLT_FLOAT_DEPTH, and BMLT_FLOAT_MISOMIP !TODO - Call separate subroutines for each of these options? call glissade_basal_melting_float(model%options%whichbmlt_float, & diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index b7f30f14..7b42087b 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/libglissade/glissade_bmlt_float.F90 @@ -41,7 +41,7 @@ module glissade_bmlt_float use glimmer_log use glide_types use cism_parallel, only: this_rank, main_task, nhalo, & - parallel_type, parallel_halo, parallel_globalindex, & + parallel_type, parallel_halo, parallel_globalindex, parallel_boundary_value, & parallel_reduce_sum, parallel_reduce_min, parallel_reduce_max implicit none @@ -52,7 +52,6 @@ module glissade_bmlt_float basin_sum, basin_average logical :: verbose_bmlt_float = .false. -!! logical :: verbose_bmlt_float = .true. logical :: verbose_velo = .true. logical :: verbose_continuity = .true. @@ -827,7 +826,9 @@ subroutine glissade_bmlt_float_thermal_forcing_init(model, ocean_data) ! This melt rate can be subtracted from the runtime melt rate to give a runtime anomaly. ! Note: On restart, bmlt_float_baseline is read from the restart file. - if (verbose_bmlt_float .and. main_task) print*, 'Compute baseline bmlt_float at initialization' + if (verbose_bmlt_float .and. main_task) then + print*, 'Compute baseline bmlt_float at initialization' + endif ! Compute some masks !TODO: Modify glissade_get_masks so that 'parallel' is not needed @@ -843,7 +844,7 @@ subroutine glissade_bmlt_float_thermal_forcing_init(model, ocean_data) call glissade_bmlt_float_thermal_forcing(& model%options%bmlt_float_thermal_forcing_param, & - model%options%ocean_data_domain, & + model%options%ocean_data_extrapolate, & parallel, & ewn, nsn, & model%numerics%dew*len0, model%numerics%dew*len0, & ! m @@ -879,6 +880,10 @@ subroutine glissade_bmlt_float_thermal_forcing_init(model, ocean_data) if (basin_number_min < 1) then + if (verbose_bmlt_float .and. main_task) then + print*, 'Extrapolate basin numbers' + endif + call basin_number_extrapolate(& ewn, nsn, & parallel, & @@ -897,7 +902,7 @@ end subroutine glissade_bmlt_float_thermal_forcing_init subroutine glissade_bmlt_float_thermal_forcing(& bmlt_float_thermal_forcing_param, & - ocean_data_domain, & + ocean_data_extrapolate, & parallel, & nx, ny, & dew, dns, & @@ -923,7 +928,7 @@ subroutine glissade_bmlt_float_thermal_forcing(& integer, intent(in) :: & bmlt_float_thermal_forcing_param, & !> melting parameterization used to derive melt rate from thermal forcing; !> current options are quadratic and ISMIP6 local, nonlocal and nonlocal_slope - ocean_data_domain !> = 0 if TF is provided on ocean domain only; = 1 if extrapolated under ice + ocean_data_extrapolate !> = 1 if TF is to be extrapolated to sub-shelf cavities, else = 0 type(parallel_type), intent(in) :: & parallel !> info for parallel communication @@ -939,6 +944,7 @@ subroutine glissade_bmlt_float_thermal_forcing(& integer, dimension(nx,ny), intent(in) :: & ice_mask, & !> = 1 where ice is present (H > 0) else = 0 marine_connection_mask !> = 1 for cells with a marine connection to the ocean, else = 0 + !> Note: marine_connection_mask includes paths through grounded marine-based cells integer, dimension(nx,ny), intent(inout) :: & ocean_mask !> = 1 for ice-free ocean, else = 0; @@ -1002,6 +1008,11 @@ subroutine glissade_bmlt_float_thermal_forcing(& tf_anomaly, & ! local version of tf_anomaly_in tf_anomaly_basin ! local version of tf_anomaly_basin_in + ! Note: This range ought to cover all regions where ice is present, but could be modified if desired. + real(dp), parameter :: & + thermal_forcing_max = 20.d0, & ! max allowed value of thermal forcing (K) + thermal_forcing_min = -5.d0 ! min allowed value of thermal forcing (K) + !TODO - Make H0_float a config parameter? real(dp), parameter :: & H0_float = 50.d0 ! thickness scale (m) for floating ice; used to reduce weights when H < H0_float @@ -1010,7 +1021,8 @@ subroutine glissade_bmlt_float_thermal_forcing(& print*, ' ' print*, 'In subroutine glissade_bmlt_float_thermal_forcing' print*, ' bmlt_float_thermal_forcing_param =', bmlt_float_thermal_forcing_param - print*, ' ocean_data_domain =', ocean_data_domain + print*, ' ocean_data_extrapolate =', ocean_data_extrapolate + print*, ' nbasin =', ocean_data%nbasin endif if (present(tf_anomaly_in)) then @@ -1031,41 +1043,11 @@ subroutine glissade_bmlt_float_thermal_forcing(& ! Make sure thermal_forcing is up to date in halo cells. call parallel_halo(ocean_data%thermal_forcing, parallel) - !WHL - Commented out the code below because this subroutine no longer uses ocean_mask to compute marine_connection_mask. - ! If CISM mis-identifies landlocked fjord cells as marine-connected, when there is no - ! marine-connected path to these cells from cells with valid data, the code will fail. - ! One possible fix is to increase the magnitude of ocean_topg_threshold in glissade_masks; - ! another is to prevent cells with non-ocean neighbors (or neighbors of neighbors) from seeding the fill. - -! if (ocean_data_domain == DATA_OCEAN_ONLY) then - ! When coupling to POP or another ocean model, TF is received at each coupling interval, - ! with unphys_val assigned to CISM cells that are outside the POP domain. - ! Some CISM cells (e.g., in fjords) may be identified as ocean (ocean_mask = 1), - ! but not have valid TF data. We do not want to count these cells as ocean cells - ! when computing marine_connection_mask, because then cells can be identified as - ! marine-connected without having a path to valid data. - ! Assume that if level k = 1 has valid data, there is valid data through the column. -! do j = 1, ny -! do i = 1, nx -! if (ocean_data%thermal_forcing(1,i,j) == unphys_val) then -! ocean_mask(i,j) = 0 -! endif -! enddo -! enddo -! endif - - if (ocean_data_domain == DATA_CISM_OCEAN_MASK) then - ! Set the thermal forcing to have unphysical values in cells where ocean_mask = 0. - ! TF will then be extrapolated into these cells at runtime, if the cells contain ice shelf cavities. - if (main_task .and. verbose_bmlt_float) print*, 'Set TF = unphys_val where CISM ocean_mask = 0' - do j = 1, ny - do i = 1, nx - if (ocean_mask(i,j) == 0) then - ocean_data%thermal_forcing(:,i,j) = unphys_val - endif - enddo - enddo - endif + ! Insert an unphysical value at the global boundary. + ! This is done to handle the case that global_bc = no_ice, + ! which puts zeroes in global boundary cells. + ! We do not want these zeroes to be interpreted as realistic thermal_forcing values. + call parallel_boundary_value(ocean_data%thermal_forcing, unphys_val, parallel) ! Set thermal_forcing_mask ! This mask identifies cells where we could have basal melting and need valid TF data. @@ -1095,17 +1077,16 @@ subroutine glissade_bmlt_float_thermal_forcing(& call parallel_halo(thermal_forcing_mask, parallel) !----------------------------------------------- - ! If thermal forcing data is provided only over the ocean domain, - ! then extrapolate the data to ice shelf cavities. - ! Note: For POP coupling, thermal forcing is received once per mass balance time step. + ! Optionally, extrapolate the ocean forcing data to ice shelf cavities. + ! Note: For coupled modeling (OCEAN_DATA_GLAD), thermal forcing is received once per mass balance time step. + ! Typically, it is computed only on the ocean grid and needs to be extrapolated to cavities. ! During the first ice sheet dynamics time step within a mass balance time step, ! several tens of iterations typically are needed to extrapolate open-ocean values ! into large shelf cavities. ! On subsequent steps, only a few iterations are needed to extrapolate into newly floating cells. !----------------------------------------------- - if (ocean_data_domain == DATA_OCEAN_ONLY .or. & - ocean_data_domain == DATA_CISM_OCEAN_MASK) then + if (ocean_data_extrapolate == OCEAN_DATA_EXTRAPOLATE_TRUE) then if (verbose_bmlt_float .and. this_rank == rtest) then print*, ' ' @@ -1197,11 +1178,11 @@ subroutine glissade_bmlt_float_thermal_forcing(& enddo endif - else ! ocean_data_domain = DATA_OCEAN_ICE; no need to extrapolate + else ! ocean data are already given everywhere; do not extrapolate if (verbose_bmlt_float .and. this_rank == rtest) then print*, ' ' - print*, 'TF to interpolate, rank, i, j =', rtest, itest, jtest + print*, 'TF to interpolate to lsrf, rank, i, j =', rtest, itest, jtest do k = kmin_diag, kmax_diag print*, ' ' print*, 'kocn =', k @@ -1214,7 +1195,7 @@ subroutine glissade_bmlt_float_thermal_forcing(& enddo endif - endif ! ocean_data_domain + endif ! ocean_data_extrapolate !----------------------------------------------- ! Interpolate the thermal forcing to the lower ice surface. @@ -1245,17 +1226,23 @@ subroutine glissade_bmlt_float_thermal_forcing(& thermal_forcing_in, & ocean_data%thermal_forcing_lsrf) - ! Bug check: Make sure there are no negative values of thermal forcing. - ! This could happen if the data set contains negative special values - ! that are not overwritten with realistic values in cavities. - !TODO - Remove this bug check if the ocean can realistically have TF < 0. + ! Bug check: Make sure there are no extreme values of thermal forcing. + ! This could happen, for example, if the input thermal forcing has special values + ! that are not overwritten with realistic values via extrapolation. + do j = 1, ny do i = 1, nx - if (ocean_data%thermal_forcing_lsrf(i,j) < 0.0d0) then + if (ocean_data%thermal_forcing_lsrf(i,j) > thermal_forcing_max) then + call parallel_globalindex(i, j, iglobal, jglobal, parallel) + write(message,*) & + 'Ocean thermal forcing error: extreme TF at i, j, lsrf, TF =', & + iglobal, jglobal, lsrf(i,j), ocean_data%thermal_forcing_lsrf(i,j) + call write_log(message, GM_FATAL) + elseif (ocean_data%thermal_forcing_lsrf(i,j) < thermal_forcing_min) then call parallel_globalindex(i, j, iglobal, jglobal, parallel) write(message,*) & - 'Ocean thermal forcing error: negative TF at level k, i, j, lsrf, TF =', & - k, iglobal, jglobal, lsrf(i,j), ocean_data%thermal_forcing_lsrf(i,j) + 'Ocean thermal forcing error: extreme TF at i, j, lsrf, TF =', & + iglobal, jglobal, lsrf(i,j), ocean_data%thermal_forcing_lsrf(i,j) call write_log(message, GM_FATAL) endif enddo @@ -1368,8 +1355,8 @@ subroutine glissade_bmlt_float_thermal_forcing(& thermal_forcing_basin, & itest, jtest, rtest) - !WHL - For diagnostics, compute the average value of deltaT_basin each basin. - ! Note: Each cell in the basin should have this average value. + ! For diagnostics, compute the average value of deltaT_basin each basin. + ! Note: Each cell in the basin should have this average value. call basin_average(& nx, ny, & @@ -1488,7 +1475,6 @@ subroutine glissade_bmlt_float_thermal_forcing(& enddo endif - !WHL - debug ! Reduce the melt rate in cells with thin floating ice, ! to reflect that these cells are only partly ice-filled. ! Note: This code gives bmlt_float = 0 in ice-free ocean cells, @@ -1570,6 +1556,7 @@ subroutine glissade_thermal_forcing_extrapolate(& integer, dimension(nx,ny), intent(in) :: & thermal_forcing_mask, & ! = 1 where thermal forcing and bmlt_float are potentially nonzero, else = 0 marine_connection_mask ! = 1 for cells with marine connection to the ocean, else = 0 + ! Note: marine_connection_mask includes paths through grounded marine-based cells real(dp), intent(in) :: & unphys_val, & ! unphysical value given to cells/levels not yet filled @@ -1603,7 +1590,7 @@ subroutine glissade_thermal_forcing_extrapolate(& max_iter, & ! max(nx,ny) * max(ewtasks, nxtasks) local_count, & ! local counter for filled values global_count, & ! global counter for filled values - global_count_save ! globalcounter for filled values from previous iteration + global_count_save ! global counter for filled values from previous iteration integer :: i, j, k, iter integer :: iglobal, jglobal @@ -1616,9 +1603,9 @@ subroutine glissade_thermal_forcing_extrapolate(& ! Note: If thermal forcing is close to but not quite equal to unphys_val (e.g., because of roundoff error), ! it is interpreted as equal to unphys_val. - real(dp), parameter :: tf_roundoff_threshold = 1.0d0 ! roundoff error threshold for thermal_forcing (deg C) + real(dp), parameter :: tf_roundoff_threshold = 1.0d0 ! roundoff error threshold for thermal_forcing (deg K) - ! Count the number of filled levels/cells in the input thermal forcing field. + ! Count the number of filled levels/cells with valid values in the input thermal_forcing. local_count = 0 do j = 1+nhalo, ny-nhalo @@ -1911,18 +1898,23 @@ subroutine glissade_thermal_forcing_extrapolate(& call parallel_halo(thermal_forcing, parallel) + ! Insert an unphysical value at the global boundary. + ! This is done to handle the case that global_bc = no_ice, + ! which puts zeroes in global boundary cells. + ! We do not want these zeroes to be interpreted as realistic thermal_forcing values. + call parallel_boundary_value(thermal_forcing, unphys_val, parallel) + ! Every several iterations, count the total number of filled cells/levels in the global domain. ! If this number has not increased since the previous iteration, then exit the loop. ! Note: Typically there are several ice sheet dynamics time steps per mass balance time step. - ! The thermal forcing field from the ocean model (ocean_data_domain = DATA_OCEAN_ONLY) is - ! updated at the start of the mass balance time step. Just after this update, the extrapolation - ! typically requires several tens of iterations to converge. But subsequent updates within this + ! The thermal forcing field from the ocean model (for OCEAN_DATA_GLAD) is updated + ! at the start of the mass balance time step. Just after this update, the extrapolation + ! typically requires 100 or more iterations to converge. But subsequent updates within this ! mass balance time step may require < 5 iterations. For this reason, we check frequently ! for convergence in the first few steps, and then less frequently as the number of iterations grows. ! Check after iterations 1, 2, and 5, then after 10, 20, etc. if (iter == 1 .or. iter == 2 .or. iter == 5 .or. mod(iter, 10) == 0) then - local_count = 0 do j = 1+nhalo, ny-nhalo do i = 1+nhalo, nx-nhalo @@ -1955,13 +1947,8 @@ subroutine glissade_thermal_forcing_extrapolate(& enddo ! max_iter - ! Make sure thermal forcing is non-negative - where (thermal_forcing /= unphys_val) - thermal_forcing = max(thermal_forcing, 0.0d0) - endwhere - ! Bug check: - ! Make sure all levels from ktop to kbot are filled in floating cells (except lakes). + ! Make sure all levels from ktop to kbot are filled in cells with thermal_forcing_mask = 1. do j = 1, ny do i = 1, nx @@ -1970,12 +1957,12 @@ subroutine glissade_thermal_forcing_extrapolate(& if (thermal_forcing(k,i,j) == unphys_val) then call parallel_globalindex(i, j, iglobal, jglobal, parallel) print*, 'i, j, ktop, kbot =', i, j, ktop(i,j), kbot(i,j) - write(message,*) & - 'Ocean data extrapolation error: did not fill level k, i, j =', k, iglobal, jglobal + write(message,*) 'Ocean data extrapolation error: unphys value in level k, i, j:', & + k, iglobal, jglobal, thermal_forcing(k,i,j) call write_log(message, GM_FATAL) endif enddo ! k - endif ! floating and not lake + endif ! thermal_forcing_mask enddo ! i enddo ! j