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/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 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 62438cbc7f..14459f7d0a 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,37 +2490,44 @@ 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]. 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 - logical :: correct_thickness + 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 ! 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 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 +2590,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 "//& @@ -2661,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 "//& @@ -2685,25 +2707,45 @@ 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 + 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 - !### 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 - 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 @@ -2720,15 +2762,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) @@ -2874,9 +2916,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) @@ -2919,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/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..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 @@ -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 @@ -700,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 @@ -711,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 @@ -762,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, & @@ -772,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 @@ -837,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 @@ -880,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 @@ -944,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 @@ -961,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 @@ -985,17 +979,17 @@ 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 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)) ) @@ -1025,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)) @@ -1069,10 +1063,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. @@ -1099,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) ) @@ -1118,10 +1112,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. @@ -1147,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) ) @@ -1254,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 @@ -1378,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 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), &