From 6e9e48e076a77042a18346bd7f1632d5ace3d528 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 19 Dec 2022 09:37:07 -0500 Subject: [PATCH 1/5] +Revise interfaces to horiz_interp_and_extract Revised the interface to the two horiz_interp_and_extract to eliminate the unused reentrant_x and tripolar_n arguments and to rename the conversion argument to scale and make it the last mandatory argument in anticipation that we might decide to make it optional, similarly to other routines like MOM_read_data. The renaming is because we use 'conversion' to indicate how the internal representation is to be rescaled for output, most prominently in register_diag_field, whereas everywhere else where we are rescaling input to the model's internal units, we use a scale argument. This makes the convention self-consistent across the MOM6 code, and should avoid some confusion. The calls to these routines were updated in 8 places in 4 modules. In 8 places the get_param calls for TRIPOLAR_N or REENTRANT_X that are no longer needed were eliminated as were the associated internal variables. This commit also includes some additions to comments near the changes that were directly tied to this commit. All answers are bitwise identical, but there are changes to a publicly visible interface; code that tries to use the old interface will not compile. --- src/framework/MOM_horizontal_regridding.F90 | 59 ++++++++++--------- .../MOM_state_initialization.F90 | 43 +++++++------- .../MOM_tracer_initialization_from_Z.F90 | 41 +++++++------ src/ocean_data_assim/MOM_oda_driver.F90 | 21 ++++--- .../vertical/MOM_ALE_sponge.F90 | 37 ++++-------- 5 files changed, 97 insertions(+), 104 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index bbb5ae0e15..7a0ace9279 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -250,14 +250,13 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, end subroutine fill_miss_2d !> Extrapolate and interpolate from a file record -subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, recnum, G, tr_z, mask_z, & - z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & +subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr_z, mask_z, & + z_in, z_edges_in, missing_value, scale, & homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol, answer_date) character(len=*), intent(in) :: filename !< Path to file containing tracer to be !! interpolated. character(len=*), intent(in) :: varnam !< Name of tracer in file. - real, intent(in) :: conversion !< Conversion factor for tracer [CU conc-1 ~> 1] integer, intent(in) :: recnum !< Record number of tracer to be read. type(ocean_grid_type), intent(inout) :: G !< Grid object real, allocatable, dimension(:,:,:), intent(out) :: tr_z @@ -271,10 +270,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real, allocatable, dimension(:), intent(out) :: z_edges_in !< Cell grid edge values for input data [Z ~> m] real, intent(out) :: missing_value !< The missing value in the returned array, scaled - !! with conversion to avoid accidentally having valid - !! values match missing values [CU ~> conc] - logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction - logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid + !! to avoid accidentally having valid values match + !! missing values [CU ~> conc] + real, intent(in) :: scale !< Scaling factor for tracer into the internal + !! units of the model [CU conc-1 ~> 1] logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data !! to produce perfectly "flat" initial conditions real, optional, intent(in) :: m_to_Z !< A conversion factor from meters to the units @@ -311,7 +310,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians] real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians] real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians] - real :: max_lat ! The maximum latitude on the input grid [degreeN] + real :: max_lat ! The maximum latitude on the input grid [degreesN] real :: pole ! The sum of tracer values at the pole [conc] real :: max_depth ! The maximum depth of the ocean [Z ~> m] real :: npole ! The number of points contributing to the pole value [nondim] @@ -330,7 +329,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, integer :: isd, ied, jsd, jed ! data domain indices integer :: id_clock_read logical :: debug=.false. - real :: I_scale ! The inverse of the conversion factor for diagnostic output [conc CU-1 ~> 1] + real :: I_scale ! The inverse of the scale factor for diagnostic output [conc CU-1 ~> 1] real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing ! iterations that determines when to stop iterating [CU ~> conc] real :: npoints ! The number of points in an average [nondim] @@ -355,10 +354,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, is_ongrid = .false. if (present(ongrid)) is_ongrid = ongrid - dtr_iter_stop = 1.0e-3*conversion + dtr_iter_stop = 1.0e-3*scale if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol - I_scale = 1.0 / conversion + I_scale = 1.0 / scale PI_180 = atan(1.0)/45. @@ -371,6 +370,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call cpu_clock_begin(id_clock_read) + ! A note by MJH copied from elsewhere suggests that this code may be using the model connectivity + ! (e.g., reentrant or tripolar) but should use the dataset's connectivity instead. + call get_var_axes_info(trim(filename), trim(varnam), axes_info) if (allocated(z_in)) deallocate(z_in) @@ -418,7 +420,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (.not. found_attr) call MOM_error(FATAL, & "error finding missing value for " // trim(varnam) // & " in file " // trim(filename) // " in hinterp_extrap") - missing_value = conversion * missing_val_in + missing_value = scale * missing_val_in call read_attribute(trim(filename), "scale_factor", scale_factor, & varname=trim(varnam), found=found_attr) @@ -471,7 +473,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, do i=is,ie if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then mask_in(i,j) = 1.0 - tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * conversion + tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * scale else tr_in(i,j) = missing_value endif @@ -511,7 +513,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, do j=1,jdp ; do i=1,id if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then mask_in(i,j) = 1.0 - tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion + tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * scale else tr_inp(i,j) = missing_value endif @@ -596,18 +598,17 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle -subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, tr_z, mask_z, & - z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & +subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, & + z_in, z_edges_in, missing_value, scale, & homogenize, spongeOngrid, m_to_Z, & answers_2018, tr_iter_tol, answer_date) integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator type(time_type), intent(in) :: Time !< A FMS time type - real, intent(in) :: conversion !< Conversion factor for tracer. type(ocean_grid_type), intent(inout) :: G !< Grid object real, allocatable, dimension(:,:,:), intent(out) :: tr_z !< Allocatable tracer array on the horizontal - !! model grid and input-file vertical levels. [CU ~> conc] + !! model grid and input-file vertical levels [CU ~> conc] real, allocatable, dimension(:,:,:), intent(out) :: mask_z !< Allocatable tracer mask array on the horizontal !! model grid and input-file vertical levels [nondim] @@ -616,10 +617,10 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t real, allocatable, dimension(:), intent(out) :: z_edges_in !< Cell grid edge values for input data [Z ~> m] real, intent(out) :: missing_value !< The missing value in the returned array, scaled - !! with conversion to avoid accidentally having valid - !! values match missing values [CU ~> conc] - logical, intent(in) :: reentrant_x !< If true, this grid is reentrant in the x-direction - logical, intent(in) :: tripolar_n !< If true, this is a northern tripolar grid + !! to avoid accidentally having valid values match + !! missing values [CU ~> conc] + real, intent(in) :: scale !< Scaling factor for tracer into the internal + !! units of the model [CU conc-1 ~> 1] logical, optional, intent(in) :: homogenize !< If present and true, horizontally homogenize data !! to produce perfectly "flat" initial conditions logical, optional, intent(in) :: spongeOngrid !< If present and true, the sponge data are on the model grid @@ -655,7 +656,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t real, dimension(:), allocatable :: lon_in ! The longitudes in the input file [degreesE] then [radians] real, dimension(:), allocatable :: lat_in ! The latitudes in the input file [degreesN] then [radians] real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole [degreesN] then [radians] - real :: max_lat ! The maximum latitude on the input grid [degreeN] + real :: max_lat ! The maximum latitude on the input grid [degreesN] real :: pole ! The sum of tracer values at the pole [conc] real :: max_depth ! The maximum depth of the ocean [Z ~> m] real :: npole ! The number of points contributing to the pole value [nondim] @@ -673,7 +674,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t logical :: debug=.false. logical :: is_ongrid integer :: ans_date ! The vintage of the expressions and order of arithmetic to use - real :: I_scale ! The inverse of the conversion factor for diagnostic output [conc CU-1 ~> 1] + real :: I_scale ! The inverse of the scale factor for diagnostic output [conc CU-1 ~> 1] real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing ! iterations that determines when to stop iterating [CU ~> conc] real :: npoints ! The number of points in an average [nondim] @@ -699,10 +700,10 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - dtr_iter_stop = 1.0e-3*conversion + dtr_iter_stop = 1.0e-3*scale if (present(tr_iter_tol)) dtr_iter_stop = tr_iter_tol - I_scale = 1.0 / conversion + I_scale = 1.0 / scale PI_180 = atan(1.0)/45. @@ -716,7 +717,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_begin(id_clock_read) call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_val_in) - missing_value = conversion*missing_val_in + missing_value = scale*missing_val_in verbosity = MOM_get_verbosity() @@ -823,7 +824,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t do j=1,jdp ; do i=1,id if (abs(tr_inp(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then mask_in(i,j) = 1.0 - tr_inp(i,j) = tr_inp(i,j) * conversion + tr_inp(i,j) = tr_inp(i,j) * scale else tr_inp(i,j) = missing_value endif @@ -909,7 +910,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t do k=1,kd do j=js,je do i=is,ie - tr_z(i,j,k) = data_in(i,j,k) * conversion + tr_z(i,j,k) = data_in(i,j,k) * scale if (ans_date >= 20190101) mask_z(i,j,k) = 1. if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. enddo diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 62438cbc7f..c5b57f3d57 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2483,8 +2483,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real :: PI_180 ! for conversion from degrees to radians [radian degree-1] real :: Hmix_default ! The default initial mixed layer depth [Z ~> m]. real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m]. - real :: missing_value_temp ! The missing value in the input temperature field - real :: missing_value_salt ! The missing value in the input salinity field + real :: missing_value_temp ! The missing value in the input temperature field [C ~> degC] + real :: missing_value_salt ! The missing value in the input salinity field [S ~> ppt] + real :: tol_temp ! The tolerance for changes in temperature during the horizontal + ! interpolation from an input dataset [C ~> degC] + real :: tol_sal ! The tolerance for changes in salinity during the horizontal + ! interpolation from an input dataset [S ~> ppt] logical :: correct_thickness real :: h_tolerance ! A parameter that controls the tolerance when adjusting the ! thickness to fit the bathymetry [Z ~> m]. @@ -2494,19 +2498,19 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just logical :: adjust_temperature = .true. ! fit t/s to target densities real :: temp_land_fill ! A temperature value to use for land points [C ~> degC] real :: salt_land_fill ! A salinity value to use for land points [C ~> degC] - logical :: reentrant_x, tripolar_n ! data arrays - real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m] - real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3] + real, dimension(:), allocatable :: z_edges_in ! Input data interface heights or depths [Z ~> m] + real, dimension(:), allocatable :: z_in ! Input data cell heights or depths [Z ~> m] + real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3] real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [C ~> degC] real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt] real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim] - real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] + real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]. - real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor - ! relative to the surface [Z ~> m]. - integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data + real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor + ! relative to the surface [Z ~> m]. + integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa]. ! Local variables for ALE remapping @@ -2569,9 +2573,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just eos => tv%eqn_of_state - reentrant_x = .false. ; call get_param(PF, mdl, "REENTRANT_X", reentrant_x, default=.true.) - tripolar_n = .false. ; call get_param(PF, mdl, "TRIPOLAR_N", tripolar_n, default=.false.) - call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE", filename, & "The name of the z-space input file used to initialize "//& "temperatures (T) and salinities (S). If T and S are not "//& @@ -2701,6 +2702,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just !### These hard-coded constants should be made into runtime parameters temp_land_fill = 0.0*US%degC_to_C salt_land_fill = 35.0*US%ppt_to_S + tol_temp = 1.0e-3*US%degC_to_C + tol_sal = 1.0e-3*US%ppt_to_S eps_z = GV%Angstrom_Z eps_rho = 1.0e-10*US%kg_m3_to_R @@ -2720,15 +2723,15 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! to the North/South Pole past the limits of the input data, they are extrapolated using the average ! value at the northernmost/southernmost latitude. - call horiz_interp_and_extrap_tracer(tfilename, potemp_var, US%degC_to_C, 1, & - G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date, & - ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%degC_to_C) + call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1, & + G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, & + scale=US%degC_to_C, homogenize=homogenize, m_to_Z=US%m_to_Z, & + answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_temp) - call horiz_interp_and_extrap_tracer(sfilename, salin_var, US%ppt_to_S, 1, & - G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date, & - ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%ppt_to_S) + call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1, & + G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, & + scale=US%ppt_to_S, homogenize=homogenize, m_to_Z=US%m_to_Z, & + answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_sal) kd = size(z_in,1) diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 04c03a5b43..7c62ea496e 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -41,11 +41,12 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. - real, dimension(:,:,:), pointer :: tr !< Pointer to array to be initialized + real, dimension(:,:,:), pointer :: tr !< Pointer to array to be initialized [CU ~> conc] type(param_file_type), intent(in) :: PF !< parameter file character(len=*), intent(in) :: src_file !< source filename character(len=*), intent(in) :: src_var_nam !< variable name in file - real, optional, intent(in) :: src_var_unit_conversion !< optional multiplicative unit conversion + real, optional, intent(in) :: src_var_unit_conversion !< optional multiplicative unit conversion, + !! often used for rescaling into model units [CU conc-1 ~> 1] integer, optional, intent(in) :: src_var_record !< record to read for multiple time-level files logical, optional, intent(in) :: homogenize !< optionally homogenize to mean value logical, optional, intent(in) :: useALEremapping !< to remap or not (optional) @@ -53,11 +54,11 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ character(len=*), optional, intent(in) :: src_var_gridspec !< Source variable name in a gridspec file. !! This is not implemented yet. ! Local variables - real :: land_fill = 0.0 - real :: convert + real :: land_fill = 0.0 ! A value to use to replace missing values [CU ~> conc] + real :: convert ! A conversion factor into the model's internal units [CU conc-1 ~> 1] integer :: recnum character(len=64) :: remapScheme - logical :: homog,useALE + logical :: homog, useALE ! This include declares and sets the variable "version". # include "version_variable.h" @@ -66,8 +67,12 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ integer :: is, ie, js, je, nz ! compute domain indices integer :: isd, ied, jsd, jed ! data domain indices integer :: i, j, k, kd - real, allocatable, dimension(:,:,:), target :: tr_z, mask_z - real, allocatable, dimension(:), target :: z_edges_in, z_in + real, allocatable, dimension(:,:,:), target :: tr_z ! Tracer array on the horizontal model grid + ! and input-file vertical levels [CU ~> conc] + real, allocatable, dimension(:,:,:), target :: mask_z ! Missing value mask on the horizontal model grid + ! and input-file vertical levels [nondim] + real, allocatable, dimension(:), target :: z_edges_in ! Cell edge depths for input data [Z ~> m] + real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m] ! Local variables for ALE remapping real, dimension(:,:,:), allocatable :: hSrc ! Source thicknesses [H ~> m or kg m-2]. @@ -75,8 +80,8 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ real :: zTopOfCell, zBottomOfCell, z_bathy ! Heights [Z ~> m]. type(remapping_CS) :: remapCS ! Remapping parameters and work arrays - real :: missing_value - integer :: nPoints + real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc] + integer :: nPoints ! The number of valid input data points in a column integer :: id_clock_routine, id_clock_ALE integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. @@ -94,7 +99,6 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ ! for horizontal regridding. Values below 20190101 recover the ! answers from 2018, while higher values use expressions that have ! been rearranged for rotational invariance. - logical :: reentrant_x, tripolar_n id_clock_routine = cpu_clock_id('(Initialize tracer from Z)', grain=CLOCK_ROUTINE) id_clock_ALE = cpu_clock_id('(Initialize tracer from Z) ALE', grain=CLOCK_LOOP) @@ -153,22 +157,17 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ "If both HOR_REGRID_2018_ANSWERS and HOR_REGRID_ANSWER_DATE are specified, the "//& "latter takes precedence.", default=default_hor_reg_ans_date) - ! These are model grid properties, but being applied to the data grid for now. - ! need to revisit this (mjh) - reentrant_x = .false. ; call get_param(PF, mdl, "REENTRANT_X", reentrant_x,default=.true.) - tripolar_n = .false. ; call get_param(PF, mdl, "TRIPOLAR_N", tripolar_n, default=.false.) - if (PRESENT(homogenize)) homog=homogenize if (PRESENT(useALEremapping)) useALE=useALEremapping if (PRESENT(remappingScheme)) remapScheme=remappingScheme - recnum=1 + recnum = 1 if (PRESENT(src_var_record)) recnum = src_var_record - convert=1.0 + convert = 1.0 if (PRESENT(src_var_unit_conversion)) convert = src_var_unit_conversion - call horiz_interp_and_extrap_tracer(src_file, src_var_nam, convert, recnum, & - G, tr_z, mask_z, z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & - homog, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date) + call horiz_interp_and_extrap_tracer(src_file, src_var_nam, recnum, & + G, tr_z, mask_z, z_in, z_edges_in, missing_value, & + scale=convert, homogenize=homog, m_to_Z=US%m_to_Z, answer_date=hor_regrid_answer_date) kd = size(z_edges_in,1)-1 call pass_var(tr_z,G%Domain) @@ -221,7 +220,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ ! Fill land values do k=1,nz ; do j=js,je ; do i=is,ie if (tr(i,j,k) == missing_value) then - tr(i,j,k)=land_fill + tr(i,j,k) = land_fill endif enddo ; enddo ; enddo diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 2a1a96168a..4ad11592f9 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -559,19 +559,22 @@ subroutine get_bias_correction_tracer(Time, US, CS) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ODA_CS), pointer :: CS !< ocean DA control structure - integer :: i,j,k + ! Local variables real, allocatable, dimension(:,:,:) :: T_bias ! Temperature biases [C ~> degC] real, allocatable, dimension(:,:,:) :: S_bias ! Salinity biases [C ~> degC] - real, allocatable, dimension(:,:,:) :: mask_z - real, allocatable, dimension(:), target :: z_in, z_edges_in - real :: missing_value - integer,dimension(3) :: fld_sz + real, allocatable, dimension(:,:,:) :: mask_z ! Missing value mask on the horizontal model grid + ! and input-file vertical levels [nondim] + real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m] + real, allocatable, dimension(:), target :: z_edges_in ! Cell edge depths for input data [Z ~> m] + real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc] + integer, dimension(3) :: fld_sz + integer :: i,j,k call cpu_clock_begin(id_clock_bias_adjustment) - call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id, Time, US%degC_to_C, CS%G, T_bias, & - mask_z, z_in, z_edges_in, missing_value, .true., .false., .false., .true.) - call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id, Time, US%ppt_to_S, CS%G, S_bias, & - mask_z, z_in, z_edges_in, missing_value, .true., .false., .false., .true.) + call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id, Time, CS%G, T_bias, & + mask_z, z_in, z_edges_in, missing_value, scale=US%degC_to_C, spongeOngrid=.true.) + call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id, Time, CS%G, S_bias, & + mask_z, z_in, z_edges_in, missing_value, scale=US%ppt_to_S, spongeOngrid=.true.) ! This should be replaced to use mask_z instead of the following lines ! which are intended to zero land values using an arbitrary limit. diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 9f5241bb9a..a121b05d30 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -133,9 +133,6 @@ module MOM_ALE_sponge logical :: time_varying_sponges !< True if using newer sponge code logical :: spongeDataOngrid !< True if the sponge data are on the model horizontal grid - logical :: reentrant_x !< grid is reentrant in the x direction - logical :: tripolar_N !< grid is folded at its north edge - !>@{ Diagnostic IDs integer, dimension(MAX_FIELDS_) :: id_sp_tendency !< Diagnostic ids for tracers !! tendency due to sponges @@ -257,11 +254,6 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "while later versions add parentheses for rotational symmetry. "//& "If both HOR_REGRID_2018_ANSWERS and HOR_REGRID_ANSWER_DATE are specified, the "//& "latter takes precedence.", default=default_hor_reg_ans_date) - call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & - "If true, the domain is zonally reentrant.", default=.true.) - call get_param(param_file, mdl, "TRIPOLAR_N", CS%tripolar_N, & - "Use tripolar connectivity at the northern edge of the "//& - "domain. With TRIPOLAR_N, NIGLOBAL must be even.", default=.false.) CS%time_varying_sponges = .false. CS%nz = GV%ke @@ -551,11 +543,6 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest "When defined, the incoming sponge data are "//& "assumed to be on the model grid " , & default=.false.) - call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & - "If true, the domain is zonally reentrant.", default=.true.) - call get_param(param_file, mdl, "TRIPOLAR_N", CS%tripolar_N, & - "Use tripolar connectivity at the northern edge of the "//& - "domain. With TRIPOLAR_N, NIGLOBAL must be even.", default=.false.) CS%time_varying_sponges = .true. CS%nz = GV%ke @@ -985,10 +972,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, CS%Ref_val(m)%scale, G, sp_val, & - mask_z, z_in, z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., & - spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & - answer_date=CS%hor_regrid_answer_date) + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, G, sp_val, & + mask_z, z_in, z_edges_in, missing_value, & + scale=CS%Ref_val(m)%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & + answer_date=CS%hor_regrid_answer_date) allocate( hsrc(nz_data) ) allocate( tmpT1d(nz_data) ) do c=1,CS%num_col @@ -1069,10 +1056,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then nz_data = CS%Ref_val_u%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, CS%Ref_val_u%scale, G, sp_val, & - mask_z, z_in, z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., & - spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& - answer_date=CS%hor_regrid_answer_date) + call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, G, sp_val, & + mask_z, z_in, z_edges_in, missing_value, & + scale=CS%Ref_val_u%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & + answer_date=CS%hor_regrid_answer_date) ! Initialize mask_z halos to zero before pass_var, in case of no update mask_z(G%isc-1, G%jsc:G%jec, :) = 0. @@ -1118,10 +1105,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) deallocate(sp_val, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, CS%Ref_val_v%scale, G, sp_val, & - mask_z, z_in, z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., & - spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& - answer_date=CS%hor_regrid_answer_date) + call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, G, sp_val, & + mask_z, z_in, z_edges_in, missing_value, & + scale=CS%Ref_val_v%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& + answer_date=CS%hor_regrid_answer_date) ! Initialize mask_z halos to zero before pass_var, in case of no update mask_z(G%isc:G%iec, G%jsc-1, :) = 0. mask_z(G%isc:G%iec, G%jec+1, :) = 0. From 76bb1276849747c42516f708e2fadb6867bd5763 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 19 Dec 2022 18:09:15 -0500 Subject: [PATCH 2/5] +Add 7 runtime variables for hard-coded tolerances Added the runtime variables DZ_BOTTOM_TOLERANCE, TRIM_IC_Z_TOLERANCE, DENSITY_INTERP_TOLERANCE, HORIZ_INTERP_TOL_TEMP, HORIZ_INTERP_TOL_SALIN, LAND_FILL_TEMP and LAND_FILL_SALIN to replace hard-coded dimensional constants in the routines MOM_temp_salt_initialize_from_Z, trim_for_ice, and initialize_thickness_from_file in MOM_state_initialization.F90. By default, all answers are bitwise identical, but there are new entries in the MOM_parameter_doc files for some configurations. --- .../MOM_state_initialization.F90 | 75 ++++++++++++++----- 1 file changed, 55 insertions(+), 20 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index c5b57f3d57..5af360d775 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -691,7 +691,9 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f ! them to units of m or correct sign conventions to positive upward [various] real :: h_tolerance ! A parameter that controls the tolerance when adjusting the ! thickness to fit the bathymetry [Z ~> m]. - integer :: inconsistent = 0 + real :: tol_dz_bot ! A tolerance for detecting inconsistent bottom depths when + ! correct_thickness is false [Z ~> m] + integer :: inconsistent ! The total number of cells with in consistent topography and layer thicknesses. logical :: correct_thickness character(len=40) :: mdl = "initialize_thickness_from_file" ! This subroutine's name. character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path @@ -738,6 +740,11 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", & units="m", default=0.1, scale=US%m_to_Z, do_not_log=just_read) endif + call get_param(param_file, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, & + "A tolerance for detecting inconsist topography and input layer "//& + "ticknesses when ADJUST_THICKNESS is false.", & + units="m", default=1.0, scale=US%m_to_Z, & + do_not_log=(just_read.or.correct_thickness)) call get_param(param_file, mdl, "INTERFACE_IC_VAR", eta_var, & "The variable name for initial conditions for interface heights "//& "relative to mean sea level, positive upward unless otherwise rescaled.", & @@ -762,8 +769,9 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f endif enddo ; enddo ; enddo + inconsistent = 0 do j=js,je ; do i=is,ie - if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > 1.0*US%m_to_Z) & + if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > tol_dz_bot) & inconsistent = inconsistent + 1 enddo ; enddo call sum_across_PEs(inconsistent) @@ -1188,6 +1196,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path real :: scale_factor ! A file-dependent scaling factor for the input pressure. real :: min_thickness ! The minimum layer thickness, recast into Z units [Z ~> m]. + real :: z_tolerance ! The tolerance with which to find the depth matching a specified pressure [Z ~> m]. integer :: i, j, k integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. @@ -1217,6 +1226,11 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) units="file dependent", default=1., do_not_log=just_read) call get_param(PF, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', & units='m', default=1.e-3, scale=US%m_to_Z, do_not_log=just_read) + call get_param(PF, mdl, "TRIM_IC_Z_TOLERANCE", z_tolerance, & + "The tolerance with which to find the depth matching the specified "//& + "surface pressure with TRIM_IC_FOR_P_SURF.", & + units="m", default=1.0e-5, scale=US%m_to_Z, do_not_log=just_read) + call get_param(PF, mdl, "TRIMMING_USES_REMAPPING", use_remapping, & 'When trimming the column, also remap T and S.', & default=.false., do_not_log=just_read) @@ -1270,7 +1284,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, & min_thickness, tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), & tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:), remap_CS, & - z_tol=1.0e-5*US%m_to_Z, remap_answer_date=remap_answer_date) + z_tol=z_tolerance, remap_answer_date=remap_answer_date) enddo ; enddo end subroutine trim_for_ice @@ -2476,7 +2490,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just integer :: i, j, k, ks integer :: nkml ! The number of layers in the mixed layer. - integer :: kd, inconsistent + integer :: inconsistent ! The total number of cells with in consistent topography and layer thicknesses. + integer :: kd ! The number of levels in the input data integer :: nkd ! number of levels to use for regridding input arrays real :: eps_Z ! A negligibly thin layer thickness [Z ~> m]. real :: eps_rho ! A negligibly small density difference [R ~> kg m-3]. @@ -2489,9 +2504,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! interpolation from an input dataset [C ~> degC] real :: tol_sal ! The tolerance for changes in salinity during the horizontal ! interpolation from an input dataset [S ~> ppt] - logical :: correct_thickness + logical :: correct_thickness ! If true, correct the column thicknesses to match the topography real :: h_tolerance ! A parameter that controls the tolerance when adjusting the ! thickness to fit the bathymetry [Z ~> m]. + real :: tol_dz_bot ! A tolerance for detecting inconsistent bottom depths when + ! correct_thickness is false [Z ~> m] character(len=40) :: potemp_var, salin_var integer, parameter :: niter=10 ! number of iterations for t/s adjustment to layer density @@ -2662,12 +2679,16 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "If true, all mass below the bottom removed if the "//& "topography is shallower than the thickness input file "//& "would indicate.", default=.false., do_not_log=just_read) - if (correct_thickness) then - call get_param(PF, mdl, "THICKNESS_TOLERANCE", h_tolerance, & + call get_param(PF, mdl, "THICKNESS_TOLERANCE", h_tolerance, & "A parameter that controls the tolerance when adjusting the "//& "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", & - units="m", default=0.1, scale=US%m_to_Z, do_not_log=just_read) - endif + units="m", default=0.1, scale=US%m_to_Z, & + do_not_log=(just_read.or..not.correct_thickness)) + call get_param(PF, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, & + "A tolerance for detecting inconsist topography and input layer "//& + "ticknesses when ADJUST_THICKNESS is false.", & + units="m", default=1.0, scale=US%m_to_Z, & + do_not_log=(just_read.or.correct_thickness)) call get_param(PF, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, & "If true, all the interior layers are adjusted to "//& @@ -2686,27 +2707,41 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "The mixed layer depth in the initial conditions when Z_INIT_SEPARATE_MIXED_LAYER "//& "is set to true.", units="m", default=US%Z_to_m*Hmix_default, scale=US%m_to_Z, & do_not_log=(just_read .or. .not.separate_mixed_layer)) + ! Reusing MINIMUM_DEPTH for the default mixed layer depth may be a strange choice, but + ! it reproduces previous answers. + call get_param(PF, mdl, "DENSITY_INTERP_TOLERANCE", eps_rho, & + "A small density tolerance used when finding depths in a density profile.", & + units="kg m-3", default=1.0e-10, scale=US%kg_m3_to_R, & + do_not_log=useALEremapping.or.just_read) call get_param(PF, mdl, "LAYER_Z_INIT_IC_EXTRAP_BUG", density_extrap_bug, & "If true use an expression with a vertical indexing bug for extrapolating the "//& "densities at the bottom of unstable profiles from data when finding the "//& "initial interface locations in layered mode from a dataset of T and S.", & default=.false., do_not_log=just_read) - ! Reusing MINIMUM_DEPTH for the default mixed layer depth may be a strange choice, but - ! it reproduces previous answers. endif + call get_param(PF, mdl, "LAND_FILL_TEMP", temp_land_fill, & + "A value to use to fill in ocean temperatures on land points.", & + units="degC", default=0.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "LAND_FILL_SALIN", salt_land_fill, & + "A value to use to fill in ocean salinities on land points.", & + units="1e-3", default=35.0, scale=US%ppt_to_S, do_not_log=just_read) + call get_param(PF, mdl, "HORIZ_INTERP_TOL_TEMP", tol_temp, & + "The tolerance in temperature changes between iterations when interpolating "//& + "ifrom an nput dataset using horiz_interp_and_extrap_tracer. This routine "//& + "converges slowly, so an overly small tolerance can get expensive.", & + units="degC", default=1.0e-3, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "HORIZ_INTERP_TOL_SALIN", tol_sal, & + "The tolerance in salinity changes between iterations when interpolating "//& + "ifrom an nput dataset using horiz_interp_and_extrap_tracer. This routine "//& + "converges slowly, so an overly small tolerance can get expensive.", & + units="1e-3", default=1.0e-3, scale=US%ppt_to_S, do_not_log=just_read) + if (just_read) then call cpu_clock_end(id_clock_routine) return ! All run-time parameters have been read, so return. endif - !### These hard-coded constants should be made into runtime parameters - temp_land_fill = 0.0*US%degC_to_C - salt_land_fill = 35.0*US%ppt_to_S - tol_temp = 1.0e-3*US%degC_to_C - tol_sal = 1.0e-3*US%ppt_to_S - eps_z = GV%Angstrom_Z - eps_rho = 1.0e-10*US%kg_m3_to_R ! Read input grid coordinates for temperature and salinity field ! in z-coordinate dataset. The file is REQUIRED to contain the @@ -2877,9 +2912,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just h(i,j,k) = GV%Z_to_H * (zi(i,j,K) - zi(i,j,K+1)) endif enddo ; enddo ; enddo - inconsistent=0 + inconsistent = 0 do j=js,je ; do i=is,ie - if (abs(zi(i,j,nz+1) - Z_bottom(i,j)) > 1.0*US%m_to_Z) & + if (abs(zi(i,j,nz+1) - Z_bottom(i,j)) > tol_dz_bot) & inconsistent = inconsistent + 1 enddo ; enddo call sum_across_PEs(inconsistent) From ef90f589d42c49408878316a0924d01184817328 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 19 Dec 2022 18:10:30 -0500 Subject: [PATCH 3/5] Document the units of variables in MOM_ALE_sponge Added to the comments describing a number of the internal variables in the MOM_ALE_sponge code, although given that much of this works on variables with arbitrary units, many of the units descriptions have to be simply [various]. All answers and output are bitwise identical. --- .../vertical/MOM_ALE_sponge.F90 | 78 +++++++++++-------- 1 file changed, 44 insertions(+), 34 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index a121b05d30..740c42f16a 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -69,8 +69,8 @@ module MOM_ALE_sponge integer :: id !< id for FMS external time interpolator integer :: nz_data !< The number of vertical levels in the input field. integer :: num_tlevs !< The number of time records contained in the file - real, dimension(:,:,:), pointer :: p => NULL() !< pointer to the data. - real, dimension(:,:,:), pointer :: h => NULL() !< pointer to the data grid. + real, dimension(:,:,:), pointer :: p => NULL() !< pointer to the data [various] + real, dimension(:,:,:), pointer :: h => NULL() !< pointer to the data grid [H ~> m or kg m-2] end type p3d !> A structure for creating arrays of pointers to 2D arrays with extra gridding information @@ -79,8 +79,8 @@ module MOM_ALE_sponge integer :: nz_data !< The number of vertical levels in the input field integer :: num_tlevs !< The number of time records contained in the file real :: scale = 1.0 !< A multiplicative factor by which to rescale input data - real, dimension(:,:), pointer :: p => NULL() !< pointer the data. - real, dimension(:,:), pointer :: h => NULL() !< pointer the data grid. + real, dimension(:,:), pointer :: p => NULL() !< pointer to the data [various] + real, dimension(:,:), pointer :: h => NULL() !< pointer the data grid [H ~> m or kg m-2] character(len=:), allocatable :: name !< The name of the input field character(len=:), allocatable :: long_name !< The long name of the input field character(len=:), allocatable :: unit !< The unit of the input field @@ -687,9 +687,9 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, & type(ALE_sponge_CS), pointer :: CS !< ALE sponge control structure (in/out). real, dimension(SZI_(G),SZJ_(G),CS%nz_data), & intent(in) :: sp_val !< Field to be used in the sponge, it can have an - !! arbitrary number of layers. + !! arbitrary number of layers [various] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - target, intent(in) :: f_ptr !< Pointer to the field to be damped + target, intent(in) :: f_ptr !< Pointer to the field to be damped [various] character(len=*), intent(in) :: sp_name !< The name of the tracer field character(len=*), optional, & intent(in) :: sp_long_name !< The long name of the tracer field @@ -698,9 +698,10 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, GV, f_ptr, CS, & intent(in) :: sp_unit !< The unit of the tracer field !! if not given, use the none real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any - !! contributions due to dimensional rescaling. The default is 1. + !! contributions due to dimensional rescaling [various ~> 1]. + !! The default is 1. - real :: scale_fac ! A factor by which to scale sp_val before storing it. + real :: scale_fac ! A factor by which to scale sp_val before storing it [various ~> 1] integer :: k, col character(len=256) :: mesg ! String for error messages character(len=256) :: long_name ! The long name of the tracer field @@ -749,7 +750,7 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - target, intent(in) :: f_ptr !< Pointer to the field to be damped (in). + target, intent(in) :: f_ptr !< Pointer to the field to be damped (in) [various]. type(ALE_sponge_CS), pointer :: CS !< Sponge control structure (in/out). character(len=*), intent(in) :: sp_name !< The name of the tracer field character(len=*), optional, & @@ -759,7 +760,8 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, intent(in) :: sp_unit !< The unit of the tracer field !! if not given, use 'none' real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any - !! contributions due to dimensional rescaling. The default is 1. + !! contributions due to dimensional rescaling [various ~> 1]. + !! The default is 1. ! Local variables integer :: isd, ied, jsd, jed @@ -824,9 +826,10 @@ subroutine set_up_ALE_sponge_vel_field_fixed(u_val, v_val, G, GV, u_ptr, v_ptr, real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any - !! contributions due to dimensional rescaling. The default is 1. + !! contributions due to dimensional rescaling [various ~> 1]. + !! The default is 1. - real :: scale_fac + real :: scale_fac ! A dimensional rescaling factor [various ~> 1] integer :: k, col if (.not.associated(CS)) return @@ -867,8 +870,9 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any - !! contributions due to dimensional rescaling. For varying - !! velocities the default is the same using US%m_s_to_L_T. + !! contributions due to dimensional rescaling, often in + !! [L s T-1 m-1 ~> 1]. For varying velocities the + !! default is the same as using US%m_s_to_L_T. ! Local variables logical :: override @@ -931,16 +935,19 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real :: damp ! The timestep times the local damping coefficient [nondim]. real :: I1pdamp ! I1pdamp is 1/(1 + damp). [nondim]. - real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid - real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid + real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid [various] + real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid [various] real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields - real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts - real, allocatable, dimension(:,:,:) :: mask_u ! A temporary array for field mask at u pts - real, allocatable, dimension(:,:,:) :: mask_v ! A temporary array for field mask at v pts - real, allocatable, dimension(:,:,:) :: tmp !< A temporary array for thermodynamic sponge tendency diagnostics, + real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields [various] + real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts [nondim] + real, allocatable, dimension(:,:,:) :: mask_u ! A temporary array for field mask at u pts [nondim] + real, allocatable, dimension(:,:,:) :: mask_v ! A temporary array for field mask at v pts [nondim] + real, allocatable, dimension(:,:,:) :: tmp !< A temporary array for thermodynamic sponge tendency + !! diagnostics [various] real, allocatable, dimension(:,:,:) :: tmp_u !< A temporary array for u sponge acceleration diagnostics + !! first in [L T-1 ~> m s-1] then in [L T-2 ~> m s-2] real, allocatable, dimension(:,:,:) :: tmp_v !< A temporary array for v sponge acceleration diagnostics + !! first in [L T-1 ~> m s-1] then in [L T-2 ~> m s-2] real, dimension(:), allocatable :: hsrc ! Source thicknesses [Z ~> m]. ! Local variables for ALE remapping real, dimension(:), allocatable :: tmpT1d @@ -948,12 +955,12 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real, allocatable, dimension(:), target :: z_in ! The depths (positive downward) in the input file [Z ~> m] real, allocatable, dimension(:), target :: z_edges_in ! The depths (positive downward) of the ! edges in the input file [Z ~> m] - real :: missing_value + real :: missing_value ! The missing value in the input data field [various] real :: Idt ! The inverse of the timestep [T-1 ~> s-1] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m]. - real :: sp_val_u ! Interpolation of sp_val to u-points - real :: sp_val_v ! Interpolation of sp_val to v-points + real :: sp_val_u ! Interpolation of sp_val to u-points, often a velocity in [L T-1 ~> m s-1] + real :: sp_val_v ! Interpolation of sp_val to v-points, often a velocity in [L T-1 ~> m s-1] integer :: nPoints is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -982,7 +989,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) i = CS%col_i(c) ; j = CS%col_j(c) CS%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data) ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 ; hsrc(:) = 0.0 ; tmpT1d(:) = -99.9 do k=1,nz_data if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(CS%col_i(c),CS%col_j(c)) ) @@ -1012,7 +1019,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) enddo endif - tmp_val1(:)=0.0;h_col(:)=0.0 + tmp_val1(:) = 0.0 ; h_col(:) = 0.0 do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data allocate(tmp_val2(CS%Ref_val(m)%nz_data)) @@ -1086,7 +1093,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) CS%Ref_val_u%p(1:nz_data,c) = 0.0 endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 ; hsrc(:) = 0.0 do k=1,nz_data if (mask_u(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(i,j) ) @@ -1134,7 +1141,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) CS%Ref_val_v%p(1:nz_data,c) = 0.0 endif ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0 + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 ; hsrc(:) = 0.0 do k=1,nz_data if (mask_v(i,j,k) == 1.0) then zBottomOfCell = -min( z_edges_in(k+1) - G%Z_ref, G%bathyT(i,j) ) @@ -1241,10 +1248,13 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, GV, turns, param_file) ! 3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:) ! All the index adjustment should follow from the Iresttime rotation - real, dimension(:,:), allocatable :: Iresttime_in, Iresttime - real, dimension(:,:,:), allocatable :: data_h_in, data_h - real, dimension(:,:,:), allocatable :: sp_val_in, sp_val - real, dimension(:,:,:), pointer :: sp_ptr => NULL() + real, dimension(:,:), allocatable :: Iresttime_in ! Restoring rate on the input sponges [T-1 ~> s-1] + real, dimension(:,:), allocatable :: Iresttime ! Restoring rate on the output sponges [T-1 ~> s-1] + real, dimension(:,:,:), allocatable :: data_h_in ! Grid for the input sponges [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: data_h ! Grid for the output sponges [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: sp_val_in ! Target data for the input sponges [various] + real, dimension(:,:,:), allocatable :: sp_val ! Target data for the output sponges [various] + real, dimension(:,:,:), pointer :: sp_ptr => NULL() ! Target data for the input sponges [various] integer :: c, c_i, c_j integer :: k, nz_data integer :: n @@ -1365,11 +1375,11 @@ end subroutine rotate_ALE_sponge subroutine update_ALE_sponge_field(sponge, p_old, G, GV, p_new) type(ALE_sponge_CS), intent(inout) :: sponge !< ALE sponge control struct real, dimension(:,:,:), & - target, intent(in) :: p_old !< The previous array of target values + target, intent(in) :: p_old !< The previous array of target values [various] type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - target, intent(in) :: p_new !< The new array of target values + target, intent(in) :: p_new !< The new array of target values [various] integer :: n From 0bda74df774f8aad16461c0b33399b9cafe7b934 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 20 Dec 2022 15:32:55 -0500 Subject: [PATCH 4/5] +Add 11 runtime params for determine_temperature Added 11 new runtime parameters (DETERMINE_TEMP_ADJUST_T_AND_S, DETERMINE_TEMP_T_MIN, DETERMINE_TEMP_T_MAX, DETERMINE_TEMP_S_MIN, DETERMINE_TEMP_S_MAX, DETERMINE_TEMP_T_TOLERANCE, DETERMINE_TEMP_S_TOLERANCE, DETERMINE_TEMP_RHO_TOLERANCE, DETERMINE_TEMP_DT_DS_WEIGHT, DETERMINE_TEMP_T_ADJ_RANGE, and DETERMINE_TEMP_S_ADJ_RANGE) to replace hard coded dimensional parameters used in the routine determine_temperature. This change also requires that a param_file_type argument and the logical argument just_read are passed to determine temperature. By default, all answers are bitwise identical but there are up to 10 new entries in the MOM_parameter_doc files for some layer-mode configurations with INIT_LAYERS_FROM_Z_FILE and FIT_TO_TARGET_DENSITY_IC set to true. --- .../MOM_state_initialization.F90 | 8 +- src/tracer/MOM_tracer_Z_init.F90 | 84 ++++++++++++++----- 2 files changed, 70 insertions(+), 22 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 5af360d775..14459f7d0a 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2737,6 +2737,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just units="1e-3", default=1.0e-3, scale=US%ppt_to_S, do_not_log=just_read) if (just_read) then + if ((.not.useALEremapping) .and. adjust_temperature) & + ! This call is just here to read and log the determine_temperature parameters + call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, 0, & + h, 0, G, GV, US, PF, just_read=.true.) call cpu_clock_end(id_clock_routine) return ! All run-time parameters have been read, so return. endif @@ -2957,8 +2961,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just if (adjust_temperature) then ! Finally adjust to target density ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1 - call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), tv%P_Ref, niter, & - h, ks, G, GV, US, eos) + call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, niter, & + h, ks, G, GV, US, PF, just_read) endif endif ! useALEremapping diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 85a858b8df..d887f5f3be 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -4,7 +4,7 @@ module MOM_tracer_Z_init ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe -! use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_io, only : MOM_read_data, get_var_sizes, read_attribute, read_variable use MOM_io, only : open_file_to_read, close_file_to_read @@ -556,8 +556,8 @@ end function find_limited_slope !> This subroutine determines the potential temperature and salinity that !! is consistent with the target density using provided initial guess -subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, GV, US, & - EOS, h_massless) +subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_start, G, GV, US, & + PF, just_read, h_massless) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -565,6 +565,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: salt !< salinity [S ~> ppt] real, dimension(SZK_(GV)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3]. + type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa]. integer, intent(in) :: niter !< maximum number of iterations integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer) @@ -572,7 +573,10 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, intent(in) :: h !< layer thickness, used only to avoid working on !! massless layers [H ~> m or kg m-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure + type(param_file_type), intent(in) :: PF !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: just_read !< If true, this call will only read + !! parameters without changing T or S. real, optional, intent(in) :: h_massless !< A threshold below which a layer is !! determined to be massless [H ~> m or kg m-2] @@ -600,29 +604,69 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, ! when old_fit is true [C ~> degC] real :: max_s_adj ! The largest permitted salinity changes with each iteration ! when old_fit is true [S ~> ppt] - logical :: adjust_salt, old_fit + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=40) :: mdl = "determine_temperature" ! This subroutine's name. + logical :: adjust_salt, fit_together integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, nz, itt is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - ! These hard coded parameters need to be set properly. - S_min = 0.5*US%ppt_to_S ; S_max = 65.0*US%ppt_to_S - T_max = 31.0*US%degC_to_C ; T_min = -2.0*US%degC_to_C - max_t_adj = 1.0*US%degC_to_C - max_s_adj = 0.5*US%ppt_to_S - tol_T = 1.0e-4*US%degC_to_C - tol_S = 1.0e-4*US%ppt_to_S - tol_rho = 1.0e-4*US%kg_m3_to_R - old_fit = .true. ! reproduces siena behavior + ! ### The algorithms of determine_temperature subroutine needs to be reexamined. - dT_dS_gauge = 10.0*US%degC_to_C*US%S_to_ppt ! 10 degC is weighted equivalently to 1 ppt. - ! ### The whole determine_temperature subroutine needs to be reexamined, both the algorithms - ! and the extensive use of hard-coded dimensional parameters. + call log_version(PF, mdl, version, "") - ! We will switch to the newer method which simultaneously adjusts + ! We should switch the default to the newer method which simultaneously adjusts ! temp and salt based on the ratio of the thermal and haline coefficients, once it is tested. + call get_param(PF, mdl, "DETERMINE_TEMP_ADJUST_T_AND_S", fit_together, & + "If true, simltaneously adjust the estimates of the temperature and salinity "//& + "based on the ratio of the thermal and haline coefficients. Otherwise try to "//& + "match the density by only adjusting temperatures within a maximum range before "//& + "revising estimates of the salinity.", default=.false., do_not_log=just_read) + ! These hard coded parameters need to be set properly. + call get_param(PF, mdl, "DETERMINE_TEMP_T_MIN", T_min, & + "The minimum temperature that can be found by determine_temperature.", & + units="degC", default=-2.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_T_MAX", T_max, & + "The maximum temperature that can be found by determine_temperature.", & + units="degC", default=31.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_S_MIN", S_min, & + "The minimum salinity that can be found by determine_temperature.", & + units="1e-3", default=0.5, scale=US%ppt_to_S, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_S_MAX", S_max, & + "The maximum salinity that can be found by determine_temperature.", & + units="1e-3", default=65.0, scale=US%ppt_to_S, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_T_TOLERANCE", tol_T, & + "The convergence tolerance for temperature in determine_temperature.", & + units="degC", default=1.0e-4, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_S_TOLERANCE", tol_S, & + "The convergence tolerance for temperature in determine_temperature.", & + units="1e-3", default=1.0e-4, scale=US%ppt_to_S, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_RHO_TOLERANCE", tol_rho, & + "The convergence tolerance for density in determine_temperature.", & + units="kg m-3", default=1.0e-4, scale=US%kg_m3_to_R, do_not_log=just_read) + if (fit_together) then + ! By default 10 degC is weighted equivalently to 1 ppt when minimizing changes. + call get_param(PF, mdl, "DETERMINE_TEMP_DT_DS_WEIGHT", dT_dS_gauge, & + "When extrapolating T & S to match the layer target densities, this "//& + "factor (in deg C / PSU) is combined with the derivatives of density "//& + "with T & S to determine what direction is orthogonal to density contours. "//& + "It could be based on a typical value of (dR/dS) / (dR/dT) in oceanic profiles.", & + units="degC PSU-1", default=10.0, scale=US%degC_to_C*US%S_to_ppt) + else + call get_param(PF, mdl, "DETERMINE_TEMP_T_ADJ_RANGE", max_t_adj, & + "The maximum amount by which the initial layer temperatures can be "//& + "modified in determine_temperature.", & + units="degC", default=1.0, scale=US%degC_to_C, do_not_log=just_read) + call get_param(PF, mdl, "DETERMINE_TEMP_S_ADJ_RANGE", max_S_adj, & + "The maximum amount by which the initial layer salinities can be "//& + "modified in determine_temperature.", & + units="1e-3", default=0.5, scale=US%ppt_to_S, do_not_log=just_read) + endif + + if (just_read) return ! All run-time parameters have been read, so return. press(:) = p_ref EOSdom(:) = EOS_domain(G%HI) @@ -643,7 +687,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, do k=k_start,nz ; do i=is,ie ! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then if (abs(rho(i,k)-R_tgt(k))>tol_rho) then - if (old_fit) then + if (.not.fit_together) then dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) else @@ -662,7 +706,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, endif enddo iter_loop - if (adjust_salt .and. old_fit) then ; do itt = 1,niter + if (adjust_salt .and. .not.fit_together) then ; do itt = 1,niter do k=1,nz call calculate_density(T(:,k), S(:,k), press, rho(:,k), EOS, EOSdom ) call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), & From 5573eb2eb1d508a2eade0a975536478fc967ab54 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 20 Dec 2022 15:36:17 -0500 Subject: [PATCH 5/5] +Add 2 runtime params for ice shelf temperatures Added the new runtime parameters INFLOW_SHELF_TEMPERATURE and MISSING_SHELF_TEMPERATURE to the ice_shelf_dynamics module to replace hard coded ice shelf temperatures. White space was also added around "=" in a number of places in this same module to align with the MOM6 style guide. By default, all answers are bitwise identical, but there are new entries in some MOM_parameter_doc files for cases that use the ice shelf code with DYNAMIC_SHELF_MASS. --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 99 ++++++++++++++---------- 1 file changed, 56 insertions(+), 43 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 2ed64359cb..552216f41d 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -147,7 +147,7 @@ module MOM_ice_shelf_dynamics logical :: moving_shelf_front !< Specify whether to advance shelf front (and calve). logical :: calve_to_mask !< If true, calve off the ice shelf when it passes the edge of a mask. real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m]. - + real :: T_shelf_missing !< An ice shelf temperature to use where there is no ice shelf [degC ~> C] real :: cg_tolerance !< The tolerance in the CG solver, relative to initial residual, that !! determines when to stop the conjugate gradient iterations [nondim]. real :: nonlinear_tolerance !< The fractional nonlinear tolerance, relative to the initial error, @@ -234,6 +234,8 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct + ! Local variables + real :: T_shelf_missing ! An ice shelf temperature to use where there is no ice shelf [degC ~> C] logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics character(len=40) :: mdl = "MOM_ice_shelf_dyn" ! This module's name. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -260,9 +262,12 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) endif if (active_shelf_dynamics) then + call get_param(param_file, mdl, "MISSING_SHELF_TEMPERATURE", T_shelf_missing, & + "An ice shelf temperature to use where there is no ice shelf.",& + units="degC", default=-10.0, scale=US%degC_to_C, do_not_log=.true.) allocate( CS%u_shelf(IsdB:IedB,JsdB:JedB), source=0.0 ) allocate( CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0 ) - allocate( CS%t_shelf(isd:ied,jsd:jed), source=-10.0*US%degC_to_C ) ! [C ~> degC] + allocate( CS%t_shelf(isd:ied,jsd:jed), source=T_shelf_missing ) ! [C ~> degC] allocate( CS%ice_visc(isd:ied,jsd:jed), source=0.0 ) allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Pa-3 s-1] allocate( CS%basal_traction(isd:ied,jsd:jed), source=0.0 ) ! [R L2 T-2 ~> Pa] @@ -329,6 +334,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! a restart file to the internal representation in this run. real :: vel_rescale ! A rescaling factor for horizontal velocities from the representation ! in a restart file to the internal representation in this run. + real :: T_shelf_bdry ! A default ice shelf temperature to use for ice flowing + ! in through open boundaries [C ~> degC] !This include declares and sets the variable "version". # include "version_variable.h" character(len=200) :: IC_file,filename,inputdir @@ -438,7 +445,13 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "if CONSTANT a constant value (for debugging).", & default="MODEL") + call get_param(param_file, mdl, "INFLOW_SHELF_TEMPERATURE", T_shelf_bdry, & + "A default ice shelf temperature to use for ice flowing in through "//& + "open boundaries.", units="degC", default=-15.0, scale=US%degC_to_C) endif + call get_param(param_file, mdl, "MISSING_SHELF_TEMPERATURE", CS%T_shelf_missing, & + "An ice shelf temperature to use where there is no ice shelf.",& + units="degC", default=-10.0, scale=US%degC_to_C) call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", CS%min_thickness_simple_calve, & "Min thickness rule for the VERY simple calving law",& units="m", default=0.0, scale=US%m_to_Z) @@ -447,7 +460,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! previously allocated for registration for restarts. if (active_shelf_dynamics) then - allocate( CS%t_bdry_val(isd:ied,jsd:jed), source=-15.0*US%degC_to_C) ! [C ~> degC] + allocate( CS%t_bdry_val(isd:ied,jsd:jed), source=T_shelf_bdry) ! [C ~> degC] allocate( CS%thickness_bdry_val(isd:ied,jsd:jed), source=0.0) allocate( CS%u_face_mask(Isdq:Iedq,Jsdq:Jedq), source=0.0) allocate( CS%v_face_mask(Isdq:Iedq,Jsdq:Jedq), source=0.0) @@ -1415,7 +1428,7 @@ subroutine ice_shelf_advect_thickness_x(CS, G, LB, time_step, hmask, h0, h_after do j=jsh,jeh ; do I=ish-1,ieh if (CS%u_face_mask(I,j) == 4.) then ! The flux itself is a specified boundary condition. uh_ice(I,j) = time_step * G%dyCu(I,j) * CS%u_flux_bdry_val(I,j) - elseif ((hmask(i,j)==1) .or. (hmask(i+1,j) == 1)) then + elseif ((hmask(i,j) == 1) .or. (hmask(i+1,j) == 1)) then u_face = 0.5 * (CS%u_shelf(I,J-1) + CS%u_shelf(I,J)) h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered. @@ -1494,7 +1507,7 @@ subroutine ice_shelf_advect_thickness_y(CS, G, LB, time_step, hmask, h0, h_after do J=jsh-1,jeh ; do i=ish,ieh if (CS%v_face_mask(i,J) == 4.) then ! The flux itself is a specified boundary condition. vh_ice(i,J) = time_step * G%dxCv(i,J) * CS%v_flux_bdry_val(i,J) - elseif ((hmask(i,j)==1) .or. (hmask(i,j+1) == 1)) then + elseif ((hmask(i,j) == 1) .or. (hmask(i,j+1) == 1)) then v_face = 0.5 * (CS%v_shelf(I-1,J) + CS%v_shelf(I,J)) h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered. @@ -1854,7 +1867,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) else - S(i,j)=ISS%h_shelf(i,j)-CS%bed_elev(i,j) + S(i,j) = ISS%h_shelf(i,j)-CS%bed_elev(i,j) endif enddo enddo @@ -2219,12 +2232,12 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, bathyT(i,j), dens_ratio, Usub, Vsub) - if (umask(I-1,J-1)==1) uret(I-1,J-1) = uret(I-1,J-1) + Usub(1,1) * basal_trac(i,j) + if (umask(I-1,J-1) == 1) uret(I-1,J-1) = uret(I-1,J-1) + Usub(1,1) * basal_trac(i,j) if (umask(I-1,J) == 1) uret(I-1,J) = uret(I-1,J) + Usub(1,2) * basal_trac(i,j) if (umask(I,J-1) == 1) uret(I,J-1) = uret(I,J-1) + Usub(2,1) * basal_trac(i,j) if (umask(I,J) == 1) uret(I,J) = uret(I,J) + Usub(2,2) * basal_trac(i,j) - if (vmask(I-1,J-1)==1) vret(I-1,J-1) = vret(I-1,J-1) + Vsub(1,1) * basal_trac(i,j) + if (vmask(I-1,J-1) == 1) vret(I-1,J-1) = vret(I-1,J-1) + Vsub(1,1) * basal_trac(i,j) if (vmask(I-1,J) == 1) vret(I-1,J) = vret(I-1,J) + Vsub(1,2) * basal_trac(i,j) if (vmask(I,J-1) == 1) vret(I,J-1) = vret(I,J-1) + Vsub(2,1) * basal_trac(i,j) if (vmask(I,J) == 1) vret(I,J) = vret(I,J) + Vsub(2,2) * basal_trac(i,j) @@ -2550,12 +2563,12 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, CS%bed_elev(i,j), & dens_ratio, Usubcontr, Vsubcontr) - if (CS%umask(I-1,J-1)==1) u_bdry_contr(I-1,J-1) = u_bdry_contr(I-1,J-1) + Usubcontr(1,1) * basal_trac(i,j) + if (CS%umask(I-1,J-1) == 1) u_bdry_contr(I-1,J-1) = u_bdry_contr(I-1,J-1) + Usubcontr(1,1) * basal_trac(i,j) if (CS%umask(I-1,J) == 1) u_bdry_contr(I-1,J) = u_bdry_contr(I-1,J) + Usubcontr(1,2) * basal_trac(i,j) if (CS%umask(I,J-1) == 1) u_bdry_contr(I,J-1) = u_bdry_contr(I,J-1) + Usubcontr(2,1) * basal_trac(i,j) if (CS%umask(I,J) == 1) u_bdry_contr(I,J) = u_bdry_contr(I,J) + Usubcontr(2,2) * basal_trac(i,j) - if (CS%vmask(I-1,J-1)==1) v_bdry_contr(I-1,J-1) = v_bdry_contr(I-1,J-1) + Vsubcontr(1,1) * basal_trac(i,j) + if (CS%vmask(I-1,J-1) == 1) v_bdry_contr(I-1,J-1) = v_bdry_contr(I-1,J-1) + Vsubcontr(1,1) * basal_trac(i,j) if (CS%vmask(I-1,J) == 1) v_bdry_contr(I-1,J) = v_bdry_contr(I-1,J) + Vsubcontr(1,2) * basal_trac(i,j) if (CS%vmask(I,J-1) == 1) v_bdry_contr(I,J-1) = v_bdry_contr(I,J-1) + Vsubcontr(2,1) * basal_trac(i,j) if (CS%vmask(I,J) == 1) v_bdry_contr(I,J) = v_bdry_contr(I,J) + Vsubcontr(2,2) * basal_trac(i,j) @@ -2610,7 +2623,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) enddo ; enddo n_g = CS%n_glen; eps_min = CS%eps_glen_min - CS%ice_visc(:,:)=1e22 + CS%ice_visc(:,:) = 1.0e22 ! Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) do j=jsc,jec ; do i=isc,iec @@ -2639,13 +2652,13 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j)) ) enddo ; enddo - if (trim(CS%ice_viscosity_compute)=="CONSTANT") then + if (trim(CS%ice_viscosity_compute) == "CONSTANT") then CS%ice_visc(i,j) = 1e15 * US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T * (G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging - elseif (trim(CS%ice_viscosity_compute)=="MODEL") then + elseif (trim(CS%ice_viscosity_compute) == "MODEL") then CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - elseif (trim(CS%ice_viscosity_compute)=="OBS") then + elseif (trim(CS%ice_viscosity_compute) == "OBS") then if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j) = CS%AGlen_visc(i,j)*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! Here CS%Aglen_visc(i,j) is the ice viscocity [Pa s-1] computed from obs and read from a file endif @@ -3008,23 +3021,23 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%u_face_mask_bdry(I-1+k,j))) case (3) - vmask(I-1+k,J-1)=3. - u_face_mask(I-1+k,j)=3. - umask(I-1+k,J)=3. - vmask(I-1+k,J)=3. - vmask(I-1+k,J)=3. + vmask(I-1+k,J-1) = 3. + u_face_mask(I-1+k,j) = 3. + umask(I-1+k,J) = 3. + vmask(I-1+k,J) = 3. + vmask(I-1+k,J) = 3. case (2) - u_face_mask(I-1+k,j)=2. + u_face_mask(I-1+k,j) = 2. case (4) - umask(I-1+k,J-1:J)=0. - vmask(I-1+k,J-1:J)=0. - u_face_mask(I-1+k,j)=4. + umask(I-1+k,J-1:J) = 0. + vmask(I-1+k,J-1:J) = 0. + u_face_mask(I-1+k,j) = 4. case (0) - umask(I-1+k,J-1:J)=0. - vmask(I-1+k,J-1:J)=0. - u_face_mask(I-1+k,j)=0. + umask(I-1+k,J-1:J) = 0. + vmask(I-1+k,J-1:J) = 0. + u_face_mask(I-1+k,j) = 0. case (1) ! stress free x-boundary - umask(I-1+k,J-1:J)=0. + umask(I-1+k,J-1:J) = 0. case default end select enddo @@ -3033,23 +3046,23 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) - vmask(I-1,J-1+k)=3. - umask(I-1,J-1+k)=3. - vmask(I,J-1+k)=3. - umask(I,J-1+k)=3. - v_face_mask(i,J-1+k)=3. + vmask(I-1,J-1+k) = 3. + umask(I-1,J-1+k) = 3. + vmask(I,J-1+k) = 3. + umask(I,J-1+k) = 3. + v_face_mask(i,J-1+k) = 3. case (2) - v_face_mask(i,J-1+k)=2. + v_face_mask(i,J-1+k) = 2. case (4) - umask(I-1:I,J-1+k)=0. - vmask(I-1:I,J-1+k)=0. - v_face_mask(i,J-1+k)=4. + umask(I-1:I,J-1+k) = 0. + vmask(I-1:I,J-1+k) = 0. + v_face_mask(i,J-1+k) = 4. case (0) - umask(I-1:I,J-1+k)=0. - vmask(I-1:I,J-1+k)=0. - v_face_mask(i,J-1+k)=0. + umask(I-1:I,J-1+k) = 0. + vmask(I-1:I,J-1+k) = 0. + v_face_mask(i,J-1+k) = 0. case (1) ! stress free y-boundary - vmask(I-1:I,J-1+k)=0. + vmask(I-1:I,J-1+k) = 0. case default end select enddo @@ -3223,7 +3236,7 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) if (ISS%h_shelf(i,j) > 0.0) then CS%t_shelf(i,j) = th_after_vflux(i,j) / ISS%h_shelf(i,j) else - CS%t_shelf(i,j) = -10.0*US%degC_to_C + CS%t_shelf(i,j) = CS%T_shelf_missing endif ! endif @@ -3234,11 +3247,11 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) else ! the ice is about to melt away in this case set thickness, area, and mask to zero ! NOTE: not mass conservative, should maybe scale salt & heat flux for this cell - CS%t_shelf(i,j) = -10.0*US%degC_to_C + CS%t_shelf(i,j) = CS%T_shelf_missing CS%tmask(i,j) = 0.0 endif elseif (ISS%hmask(i,j) == 0) then - CS%t_shelf(i,j) = -10.0*US%degC_to_C + CS%t_shelf(i,j) = CS%T_shelf_missing elseif ((ISS%hmask(i,j) == 3) .or. (ISS%hmask(i,j) == -2)) then CS%t_shelf(i,j) = CS%t_bdry_val(i,j) endif