diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index a04ee426e6..7075fb7c10 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -9,10 +9,8 @@ module MOM_surface_forcing_gfdl use MOM_constants, only : hlv, hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT -use MOM_diag_mediator, only : diag_ctrl -use MOM_diag_mediator, only : safe_alloc_ptr, time_type +use MOM_diag_mediator, only : diag_ctrl, safe_alloc_ptr, time_type use MOM_domains, only : pass_vector, pass_var, fill_symmetric_edges -use MOM_domains, only : global_field_sum, BITWISE_EXACT_SUM use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE, To_All use MOM_domains, only : To_North, To_East, Omit_Corners use MOM_error_handler, only : MOM_error, WARNING, FATAL, is_root_pe, MOM_mesg diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9b650f8598..37ebeda1fa 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -62,6 +62,7 @@ module MOM_open_boundary public update_OBC_ramp public rotate_OBC_config public rotate_OBC_init +public initialize_segment_data integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary @@ -268,7 +269,7 @@ module MOM_open_boundary real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of !! characteristics) in units of grid points per timestep [nondim]. logical :: OBC_pe !< Is there an open boundary on this tile? - type(remapping_CS), pointer :: remap_CS !< ALE remapping control structure for segments only + type(remapping_CS), pointer :: remap_CS=> NULL() !< ALE remapping control structure for segments only type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries real, pointer, dimension(:,:,:) :: & rx_normal => NULL(), & !< Array storage for normal phase speed for EW radiation OBCs in units of @@ -341,6 +342,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str" character(len=200) :: config1 ! String for OBC_USER_CONFIG real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] + character(len=128) :: inputdir + logical :: answers_2018, default_2018_answers + logical :: check_reconstruction, check_remapping, force_bounds_in_subcell + character(len=32) :: remappingScheme + allocate(OBC) call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", OBC%number_of_segments, & @@ -497,7 +503,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) enddo ! if (open_boundary_query(OBC, needs_ext_seg_data=.true.)) & - call initialize_segment_data(G, OBC, param_file) + ! call initialize_segment_data(G, OBC, param_file) if (open_boundary_query(OBC, apply_open_OBC=.true.)) then call get_param(param_file, mdl, "OBC_RADIATION_MAX", OBC%rx_max, & @@ -540,6 +546,39 @@ subroutine open_boundary_config(G, US, param_file, OBC) if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out enddo + call get_param(param_file, mdl, "REMAPPING_SCHEME", remappingScheme, & + "This sets the reconstruction scheme used "//& + "for vertical remapping for all variables. "//& + "It can be one of the following schemes: \n"//& + trim(remappingSchemesDoc), default=remappingDefaultScheme,do_not_log=.true.) + call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & + "If true, cell-by-cell reconstructions are checked for "//& + "consistency and if non-monotonicity or an inconsistency is "//& + "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, & + "If true, the results of remapping are checked for "//& + "conservation and new extrema and if an inconsistency is "//& + "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, & + "If true, read external OBC data on the supergrid.", & + default=.false.) + call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & + "If true, the values on the intermediate grid used for remapping "//& + "are forced to be bounded, which might not be the case due to "//& + "round off.", default=.false.,do_not_log=.true.) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.false.) + call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & + "If true, use the order of arithmetic and expressions that recover the "//& + "answers from the end of 2018. Otherwise, use updated and more robust "//& + "forms of the same expressions.", default=default_2018_answers) + + allocate(OBC%remap_CS) + call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & + check_reconstruction=check_reconstruction, check_remapping=check_remapping, & + force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) + endif ! OBC%number_of_segments > 0 ! Safety check @@ -564,7 +603,7 @@ end subroutine open_boundary_config subroutine initialize_segment_data(G, OBC, PF) use mpp_mod, only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes - type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure type(param_file_type), intent(in) :: PF !< Parameter file handle @@ -576,10 +615,7 @@ subroutine initialize_segment_data(G, OBC, PF) character(len=32), dimension(MAX_OBC_FIELDS) :: fields ! segment field names character(len=128) :: inputdir type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list - character(len=32) :: remappingScheme character(len=256) :: mesg ! Message for error messages. - logical :: check_reconstruction, check_remapping, force_bounds_in_subcell - logical :: answers_2018, default_2018_answers integer, dimension(4) :: siz,siz2 integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -599,39 +635,6 @@ subroutine initialize_segment_data(G, OBC, PF) call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) - call get_param(PF, mdl, "REMAPPING_SCHEME", remappingScheme, & - "This sets the reconstruction scheme used "//& - "for vertical remapping for all variables. "//& - "It can be one of the following schemes: \n"//& - trim(remappingSchemesDoc), default=remappingDefaultScheme,do_not_log=.true.) - call get_param(PF, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, & - "If true, cell-by-cell reconstructions are checked for "//& - "consistency and if non-monotonicity or an inconsistency is "//& - "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "FATAL_CHECK_REMAPPING", check_remapping, & - "If true, the results of remapping are checked for "//& - "conservation and new extrema and if an inconsistency is "//& - "detected then a FATAL error is issued.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, & - "If true, the values on the intermediate grid used for remapping "//& - "are forced to be bounded, which might not be the case due to "//& - "round off.", default=.false.,do_not_log=.true.) - call get_param(PF, mdl, "BRUSHCUTTER_MODE", OBC%brushcutter_mode, & - "If true, read external OBC data on the supergrid.", & - default=.false.) - call get_param(PF, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & - "This sets the default value for the various _2018_ANSWERS parameters.", & - default=.false.) - call get_param(PF, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & - "If true, use the order of arithmetic and expressions that recover the "//& - "answers from the end of 2018. Otherwise, use updated and more robust "//& - "forms of the same expressions.", default=default_2018_answers) - - allocate(OBC%remap_CS) - call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & - check_reconstruction=check_reconstruction, check_remapping=check_remapping, & - force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) - if (OBC%user_BCs_set_globally) return ! Try this here just for the documentation. It is repeated below. @@ -4966,6 +4969,8 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) integer :: l + if (OBC_in%number_of_segments==0) return + ! Scalar and logical transfer OBC%number_of_segments = OBC_in%number_of_segments OBC%ke = OBC_in%ke @@ -5023,8 +5028,10 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) OBC%OBC_pe = OBC_in%OBC_pe ! remap_CS is set up by initialize_segment_data, so we copy the fields here. - allocate(OBC%remap_CS) - OBC%remap_CS = OBC_in%remap_CS + if (ASSOCIATED(OBC_in%remap_CS)) then + allocate(OBC%remap_CS) + OBC%remap_CS = OBC_in%remap_CS + endif ! TODO: The OBC registry seems to be a list of "registered" OBC types. ! It does not appear to be used, so for now we skip this record. diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index e3a5443840..47a2bf21b0 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -473,7 +473,7 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom - intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) + intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) do m=2,4 wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6b68cb3deb..66fd873f67 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -51,7 +51,6 @@ module MOM_ice_shelf use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init -use time_manager_mod, only : print_time implicit none ; private #include diff --git a/src/ice_shelf/MOM_ice_shelf_state.F90 b/src/ice_shelf/MOM_ice_shelf_state.F90 index b3e88697f2..a3784b5a34 100644 --- a/src/ice_shelf/MOM_ice_shelf_state.F90 +++ b/src/ice_shelf/MOM_ice_shelf_state.F90 @@ -12,7 +12,6 @@ module MOM_ice_shelf_state use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type use MOM_get_input, only : directories, Get_MOM_input -use mpp_mod, only : mpp_sum, mpp_max, mpp_min, mpp_pe, mpp_npes, mpp_sync use MOM_coms, only : reproducing_sum use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index e451966364..b613648c7c 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -18,13 +18,11 @@ module MOM_state_initialization use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type, isPointInCell use MOM_interface_heights, only : find_eta -use MOM_io, only : file_exists -use MOM_io, only : MOM_read_data, MOM_read_vector -use MOM_io, only : slasher -use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init +use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher +use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init, set_tracer_data use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE use MOM_open_boundary, only : open_boundary_query -use MOM_open_boundary, only : set_tracer_data +use MOM_open_boundary, only : set_tracer_data, initialize_segment_data use MOM_open_boundary, only : open_boundary_test_extern_h use MOM_open_boundary, only : fill_temp_salt_segments use MOM_open_boundary, only : update_OBC_segment_data @@ -33,8 +31,7 @@ module MOM_state_initialization use MOM_restart, only : restore_state, determine_is_new_run, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density use MOM_sponge, only : initialize_sponge, sponge_CS -use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge -use MOM_ALE_sponge, only : ALE_sponge_CS +use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge, ALE_sponge_CS use MOM_string_functions, only : uppercase, lowercase use MOM_time_manager, only : time_type use MOM_tracer_registry, only : tracer_registry_type @@ -44,8 +41,7 @@ module MOM_state_initialization use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type, EOS_domain use MOM_EOS, only : convert_temp_salt_for_TEOS10 use user_initialization, only : user_initialize_thickness, user_initialize_velocity -use user_initialization, only : user_init_temperature_salinity -use user_initialization, only : user_set_OBC_data +use user_initialization, only : user_init_temperature_salinity, user_set_OBC_data use user_initialization, only : user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness use DOME_initialization, only : DOME_set_OBC_data @@ -97,7 +93,6 @@ module MOM_state_initialization use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer -use fms_io_mod, only : field_size implicit none ; private @@ -563,6 +558,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! This controls user code for setting open boundary data if (associated(OBC)) then + call initialize_segment_data(G, OBC, PF) ! call initialize_segment_data(G, OBC, param_file) +! call open_boundary_config(G, US, PF, OBC) ! Call this once to fill boundary arrays from fixed values if (.not. OBC%needs_IO_for_data) & call update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 91085047c9..713282cdb3 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -195,14 +195,14 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif hc = (tv%C_p*GV%H_to_RZ) * h(i,j,k) - if (h(i,j,k) <= 10.0*GV%Angstrom_H) then + if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then ! Very thin layers should not be cooled by the frazil flux. if (tv%T(i,j,k) < T_freeze(i)) then fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) tv%T(i,j,k) = T_freeze(i) endif - else - if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) <= 0.0) then + elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then + if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc fraz_col(i) = 0.0 else @@ -822,19 +822,18 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - ! Only apply forcing if fluxes%sw is associated. - if (.not.associated(fluxes%sw)) return - -#define _OLD_ALG_ Idt = 1.0 / dt calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) calculate_buoyancy = present(SkinBuoyFlux) if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0 + if (present(cTKE)) cTKE(:,:,:) = 0.0 g_Hconv2 = (US%L_to_Z**2*GV%g_Earth * GV%H_to_RZ) * GV%H_to_RZ EOSdom(:) = EOS_domain(G%HI) - if (present(cTKE)) cTKE(:,:,:) = 0.0 + ! Only apply forcing if fluxes%sw is associated. + if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return + if (calculate_buoyancy) then SurfPressure(:) = 0.0 GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 @@ -874,7 +873,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t h2d(i,k) = h(i,j,k) T2d(i,k) = tv%T(i,j,k) enddo ; enddo - if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) if (calculate_energetics) then ! The partial derivatives of specific volume with temperature and @@ -898,6 +896,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t pen_TKE_2d(:,:) = 0.0 endif + ! Nothing more is done on this j-slice if there is no buoyancy forcing. + if (.not.associated(fluxes%sw)) cycle + + if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) + ! The surface forcing is contained in the fluxes type. ! We aggregate the thermodynamic forcing for a time step into the following: ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c6a6f37b39..1a4fb58e02 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1143,12 +1143,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: z2 ! A copy of z_i [nondim] + real :: botfn ! A function that is 1 at the bottom and small far from it [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] logical :: do_shelf, do_OBCs integer :: i, k, is, ie, max_nk integer :: nz - real :: botfn a_cpl(:,:) = 0.0 Kv_tot(:,:) = 0.0 diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 7d2310b42f..66c0e33bac 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -13,10 +13,8 @@ module MOM_generic_tracer #define _ALLOCATED allocated #endif - ! ### These imports should not reach into FMS directly ### - use mpp_mod, only: stdout, mpp_error, FATAL,WARNING,NOTE - use field_manager_mod, only: fm_get_index,fm_string_len + use field_manager_mod, only: fm_string_len use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag @@ -108,7 +106,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Local variables logical :: register_MOM_generic_tracer - character(len=fm_string_len), parameter :: sub_name = 'register_MOM_generic_tracer' + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' character(len=200) :: inputdir ! The directory where NetCDF input files are. ! These can be overridden later in via the field manager? @@ -122,7 +120,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) register_MOM_generic_tracer = .false. if (associated(CS)) then - call mpp_error(WARNING, "register_MOM_generic_tracer called with an "// & + call MOM_error(WARNING, "register_MOM_generic_tracer called with an "// & "associated control structure.") return endif @@ -185,7 +183,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) !Get the tracer list call generic_tracer_get_list(CS%g_tracer_list) - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ": No tracer in the list.") ! For each tracer name get its T_prog index and get its fields @@ -247,7 +245,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the !! ALE sponges. - character(len=fm_string_len), parameter :: sub_name = 'initialize_MOM_generic_tracer' + character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer' logical :: OK integer :: i, j, k, isc, iec, jsc, jec, nk type(g_tracer_type), pointer :: g_tracer,g_tracer_next @@ -265,7 +263,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, CS%diag=>diag !Get the tracer list - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ": No tracer in the list.") !For each tracer name get its fields g_tracer=>CS%g_tracer_list @@ -426,7 +424,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) ! Local variables - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_column_physics' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics' type(g_tracer_type), pointer :: g_tracer, g_tracer_next character(len=fm_string_len) :: g_tracer_name @@ -443,7 +441,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = G%ke !Get the tracer list - if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL,& + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,& trim(sub_name)//": No tracer in the list.") #ifdef _USE_MOM6_DIAG @@ -587,7 +585,7 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde type(g_tracer_type), pointer :: g_tracer, g_tracer_next real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_stock' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock' integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -660,7 +658,7 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, yg type(g_tracer_type), pointer :: g_tracer, g_tracer_next real, dimension(:,:,:,:), pointer :: tr_field real, dimension(:,:,:), pointer :: tr_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_min_max' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max' real, dimension(:,:,:),pointer :: grid_tmask integer :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau @@ -728,7 +726,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, CS) ! Local variables real :: sosga - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_surface_state' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state' real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke,1) :: rho0 real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke) :: dzt type(g_tracer_type), pointer :: g_tracer @@ -750,7 +748,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, CS) tau=1,sosga=sosga,model_time=get_diag_time_end(CS%diag)) !Output diagnostics via diag_manager for all tracers in this module -! if (.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& +! if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& ! "No tracer in the list.") ! call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld @@ -767,7 +765,7 @@ subroutine MOM_generic_flux_init(verbosity) integer :: ind character(len=fm_string_len) :: g_tracer_name,longname, package,units,old_package,file_in,file_out real :: const_init_value - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_flux_init' + character(len=128), parameter :: sub_name = 'MOM_generic_flux_init' type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next if (.not. g_registered) then @@ -777,7 +775,7 @@ subroutine MOM_generic_flux_init(verbosity) call generic_tracer_get_list(g_tracer_list) if (.NOT. associated(g_tracer_list)) then - call mpp_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") + call MOM_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") return endif @@ -812,7 +810,7 @@ subroutine MOM_generic_tracer_get(name,member,array, CS) type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. real, dimension(:,:,:), pointer :: array_ptr - character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_get' + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr) array(:,:,:) = array_ptr(:,:,:) diff --git a/src/tracer/MOM_offline_aux.F90 b/src/tracer/MOM_offline_aux.F90 index 21db2cfff4..119ad555da 100644 --- a/src/tracer/MOM_offline_aux.F90 +++ b/src/tracer/MOM_offline_aux.F90 @@ -4,7 +4,6 @@ module MOM_offline_aux ! This file is part of MOM6. See LICENSE.md for the license. -use mpp_domains_mod, only : CENTER, CORNER, NORTH, EAST use data_override_mod, only : data_override_init, data_override use MOM_time_manager, only : time_type, operator(-) use MOM_debugging, only : check_column_integrals @@ -12,7 +11,7 @@ module MOM_offline_aux use MOM_diag_vkernels, only : reintegrate_column use MOM_error_handler, only : callTree_enter, callTree_leave, MOM_error, FATAL, WARNING, is_root_pe use MOM_grid, only : ocean_grid_type -use MOM_io, only : MOM_read_data, MOM_read_vector +use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_verticalGrid, only : verticalGrid_type use MOM_file_parser, only : get_param, log_version, param_file_type use astronomy_mod, only : orbital_time, diurnal_solar, daily_mean_solar diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index b7af9849b3..3895e8a116 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -4,7 +4,6 @@ module MOM_offline_main ! This file is part of MOM6. See LICENSE.md for the license. -use mpp_domains_mod, only : CENTER, CORNER, NORTH, EAST use MOM_ALE, only : ALE_CS, ALE_main_offline, ALE_offline_inputs use MOM_checksums, only : hchksum, uvchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end @@ -20,7 +19,7 @@ module MOM_offline_main use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type -use MOM_io, only : MOM_read_data, MOM_read_vector +use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_offline_aux, only : update_offline_from_arrays, update_offline_from_files use MOM_offline_aux, only : next_modulo_time, offline_add_diurnal_sw use MOM_offline_aux, only : update_h_horizontal_flux, update_h_vertical_flux, limit_mass_flux_3d diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 6a362d4fd5..e9c8fb0e7b 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -140,16 +140,15 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & !$OMP hprev,domore_k,js,je,is,ie,uhtr,vhtr,G,GV,h_end,& !$OMP uh_neglect,vh_neglect,ntr,Tr,h_prev_opt) -! This initializes the halos of uhr and vhr because pass_vector might do -! calculations on them, even though they are never used. -!$OMP do - + ! This initializes the halos of uhr and vhr because pass_vector might do + ! calculations on them, even though they are never used. + !$OMP do do k=1,nz do j=jsd,jed ; do I=IsdB,IedB ; uhr(I,j,k) = 0.0 ; enddo ; enddo do J=jsdB,jedB ; do i=Isd,Ied ; vhr(i,J,k) = 0.0 ; enddo ; enddo do j=jsd,jed ; do i=Isd,Ied ; hprev(i,j,k) = 0.0 ; enddo ; enddo domore_k(k)=1 -! Put the remaining (total) thickness fluxes into uhr and vhr. + ! Put the remaining (total) thickness fluxes into uhr and vhr. do j=js,je ; do I=is-1,ie ; uhr(I,j,k) = uhtr(I,j,k) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; vhr(i,J,k) = vhtr(i,J,k) ; enddo ; enddo if (.not. present(h_prev_opt)) then @@ -173,17 +172,17 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & enddo -!$OMP do + !$OMP do do j=jsd,jed ; do I=isd,ied-1 - uh_neglect(I,j) = GV%H_subroundoff*MIN(G%areaT(i,j),G%areaT(i+1,j)) + uh_neglect(I,j) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i+1,j)) enddo ; enddo -!$OMP do + !$OMP do do J=jsd,jed-1 ; do i=isd,ied - vh_neglect(i,J) = GV%H_subroundoff*MIN(G%areaT(i,j),G%areaT(i,j+1)) + vh_neglect(i,J) = GV%H_subroundoff * MIN(G%areaT(i,j), G%areaT(i,j+1)) enddo ; enddo -!$OMP do ! initialize diagnostic fluxes and tendencies + !$OMP do do m=1,ntr if (associated(Tr(m)%ad_x)) then do k=1,nz ; do j=jsd,jed ; do i=isd,ied @@ -207,7 +206,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & do J=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_y(i,J) = 0.0 ; enddo ; enddo endif enddo -!$OMP end parallel + !$OMP end parallel isv = is ; iev = ie ; jsv = js ; jev = je @@ -222,8 +221,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & ! Reevaluate domore_u & domore_v unless the valid range is the same size as ! before. Also, do this if there is Strang splitting. if ((nsten_halo > 1) .or. (itt==1)) then -!$OMP parallel do default(none) shared(nz,domore_k,jsv,jev,domore_u,isv,iev,stencil, & -!$OMP uhr,domore_v,vhr) + !$OMP parallel do default(shared) do k=1,nz ; if (domore_k(k) > 0) then do j=jsv,jev ; if (.not.domore_u(j,k)) then do i=isv+stencil-1,iev-stencil; if (uhr(I,j,k) /= 0.0) then @@ -256,9 +254,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & ! for all the transport to happen. The sum over domore_k keeps the processors ! synchronized. This may not be very efficient, but it should be reliable. -!$OMP parallel default(private) shared(nz,domore_k,x_first,Tr,hprev,uhr,uh_neglect, & -!$OMP OBC,domore_u,ntr,Idt,isv,iev,jsv,jev,stencil, & -!$OMP G,GV,CS,vhr,vh_neglect,domore_v,US) + !$OMP parallel default(shared) if (x_first) then @@ -305,7 +301,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, & endif ! x_first -!$OMP end parallel + !$OMP end parallel ! If the advection just isn't finishing after max_iter, move on. if (itt >= max_iter) then @@ -385,6 +381,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. real :: min_h ! The minimum thickness that can be realized during ! any of the passes [H ~> m or kg m-2]. + real :: tiny_h ! The smallest numerically invertable thickness [H ~> m or kg m-2]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points. @@ -406,16 +403,15 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & if (usePPM .and. .not. useHuynh) stencil = 2 min_h = 0.1*GV%Angstrom_H + tiny_h = tiny(min_h) h_neglect = GV%H_subroundoff -! do I=is-1,ie ; ts2(I) = 0.0 ; enddo do I=is-1,ie ; CFL(I) = 0.0 ; enddo do j=js,je ; if (domore_u(j,k)) then domore_u(j,k) = .false. - ! Calculate the i-direction profiles (slopes) of each tracer that - ! is being advected. + ! Calculate the i-direction profiles (slopes) of each tracer that is being advected. if (usePLMslope) then do m=1,ntr ; do i=is-stencil,ie+stencil !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < & @@ -490,33 +486,33 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! in the cell plus whatever part of its half of the mass flux that ! the flux through the other side does not require. do I=is-1,ie - if (uhr(I,j,k) == 0.0) then + if ((uhr(I,j,k) == 0.0) .or. & + ((uhr(I,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. & + ((uhr(I,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then uhh(I) = 0.0 CFL(I) = 0.0 elseif (uhr(I,j,k) < 0.0) then hup = hprev(i+1,j,k) - G%areaT(i+1,j)*min_h - hlos = MAX(0.0,uhr(I+1,j,k)) + hlos = MAX(0.0, uhr(I+1,j,k)) if ((((hup - hlos) + uhr(I,j,k)) < 0.0) .and. & ((0.5*hup + uhr(I,j,k)) < 0.0)) then - uhh(I) = MIN(-0.5*hup,-hup+hlos,0.0) + uhh(I) = MIN(-0.5*hup, -hup+hlos, 0.0) domore_u(j,k) = .true. else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j))) - CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive + CFL(I) = - uhh(I) / (hprev(i+1,j,k)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h - hlos = MAX(0.0,-uhr(I-1,j,k)) + hlos = MAX(0.0, -uhr(I-1,j,k)) if ((((hup - hlos) - uhr(I,j,k)) < 0.0) .and. & ((0.5*hup - uhr(I,j,k)) < 0.0)) then - uhh(I) = MAX(0.5*hup,hup-hlos,0.0) + uhh(I) = MAX(0.5*hup, hup-hlos, 0.0) domore_u(j,k) = .true. else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) - CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive + CFL(I) = uhh(I) / (hprev(i,j,k)) ! CFL is positive endif enddo @@ -545,11 +541,11 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & dA = aR - aL ; mA = 0.5*( aR + aL ) if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then - aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells + aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then - aL = 3.*Tc - 2.*aR + aL = 3.*Tc - 2.*aR elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then - aR = 3.*Tc - 2.*aL + aR = 3.*Tc - 2.*aL endif a6 = 6.*Tc - 3. * (aR + aL) ! Curvature @@ -570,28 +566,17 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m) - !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) ) - ! Alternative implementation of PLM Tc = T_tmp(i,m) flux_x(I,j,m) = uhh(I)*( Tc + 0.5 * slope_x(i,m) * ( 1. - CFL(I) ) ) - ! Original implementation of PLM - !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I)) else ! Indirect implementation of PLM !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m) !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) ) ! Alternative implementation of PLM - !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m) - !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) ) - ! Alternative implementation of PLM Tc = T_tmp(i+1,m) flux_x(I,j,m) = uhh(I)*( Tc - 0.5 * slope_x(i+1,m) * ( 1. - CFL(I) ) ) - ! Original implementation of PLM - !flux_x(I,j,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j))) enddo ; enddo endif ! usePPM @@ -760,6 +745,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & CFL ! The absolute value of the advective upwind-cell CFL number [nondim]. real :: min_h ! The minimum thickness that can be realized during ! any of the passes [H ~> m or kg m-2]. + real :: tiny_h ! The smallest numerically invertable thickness [H ~> m or kg m-2]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. @@ -777,8 +763,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & if (usePPM .and. .not. useHuynh) stencil = 2 min_h = 0.1*GV%Angstrom_H + tiny_h = tiny(min_h) h_neglect = GV%H_subroundoff - !do i=is,ie ; ts2(i) = 0.0 ; enddo ! We conditionally perform work on tracer points: calculating the PLM slope, ! and updating tracer concentration within a cell @@ -822,7 +808,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! make a copy of the tracers in case values need to be overridden for OBCs - do j=G%jsd,G%jed; do m=1,ntr; do i=G%isd,G%ied + do j=G%jsd,G%jed ; do m=1,ntr ; do i=G%isd,G%ied T_tmp(i,m,j) = Tr(m)%t(i,j,k) enddo ; enddo ; enddo @@ -873,33 +859,33 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & domore_v(J,k) = .false. do i=is,ie - if (vhr(i,J,k) == 0.0) then + if ((vhr(i,J,k) == 0.0) .or. & + ((vhr(i,J,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. & + ((vhr(i,J,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then vhh(i,J) = 0.0 CFL(i) = 0.0 elseif (vhr(i,J,k) < 0.0) then hup = hprev(i,j+1,k) - G%areaT(i,j+1)*min_h - hlos = MAX(0.0,vhr(i,J+1,k)) + hlos = MAX(0.0, vhr(i,J+1,k)) if ((((hup - hlos) + vhr(i,J,k)) < 0.0) .and. & ((0.5*hup + vhr(i,J,k)) < 0.0)) then - vhh(i,J) = MIN(-0.5*hup,-hup+hlos,0.0) + vhh(i,J) = MIN(-0.5*hup, -hup+hlos, 0.0) domore_v(J,k) = .true. else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) - CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive + CFL(i) = - vhh(i,J) / hprev(i,j+1,k) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h - hlos = MAX(0.0,-vhr(i,J-1,k)) + hlos = MAX(0.0, -vhr(i,J-1,k)) if ((((hup - hlos) - vhr(i,J,k)) < 0.0) .and. & ((0.5*hup - vhr(i,J,k)) < 0.0)) then - vhh(i,J) = MAX(0.5*hup,hup-hlos,0.0) + vhh(i,J) = MAX(0.5*hup, hup-hlos, 0.0) domore_v(J,k) = .true. else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) - CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive + CFL(i) = vhh(i,J) / hprev(i,j,k) ! CFL is positive endif enddo @@ -952,26 +938,16 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) ) ! Alternative implementation of PLM - !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j) - !flux_y(i,m,J) = vhh(i,J)*(aR - 0.5 * slope_y(i,m,j)*CFL(i)) - ! Alternative implementation of PLM Tc = T_tmp(i,m,j) flux_y(i,m,J) = vhh(i,J)*( Tc + 0.5 * slope_y(i,m,j) * ( 1. - CFL(i) ) ) - ! Original implementation of PLM - !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j,k) + slope_y(i,m,j)*ts2(i)) else ! Indirect implementation of PLM !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1) !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) ) ! Alternative implementation of PLM - !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1) - !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * slope_y(i,m,j+1)*CFL(i) ) - ! Alternative implementation of PLM Tc = T_tmp(i,m,j+1) flux_y(i,m,J) = vhh(i,J)*( Tc - 0.5 * slope_y(i,m,j+1) * ( 1. - CFL(i) ) ) - ! Original implementation of PLM - !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j+1,k) - slope_y(i,m,j+1)*ts2(i)) endif enddo ; enddo endif ! usePPM diff --git a/src/tracer/RGC_tracer.F90 b/src/tracer/RGC_tracer.F90 index 028718f379..44c6c2e5a1 100644 --- a/src/tracer/RGC_tracer.F90 +++ b/src/tracer/RGC_tracer.F90 @@ -19,7 +19,7 @@ module RGC_tracer use MOM_forcing_type, only : forcing use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc use MOM_restart, only : MOM_restart_CS use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS, get_ALE_sponge_nz_data use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -207,8 +207,7 @@ subroutine initialize_RGC_tracer(restart, day, G, GV, h, diag, OBC, CS, & CS%tracer_IC_file) do m=1,NTR call query_vardesc(CS%tr_desc(m), name, caller="initialize_RGC_tracer") - call read_data(CS%tracer_IC_file, trim(name), & - CS%tr(:,:,:,m), domain=G%Domain%mpp_domain) + call MOM_read_data(CS%tracer_IC_file, trim(name), CS%tr(:,:,:,m), G%Domain) enddo else do m=1,NTR