From f54338ad1493a25d304ff48943d234843bc3f48c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 14 Apr 2022 15:45:52 -0400 Subject: [PATCH] Fix non-standard white space (#109) * Fix non-standard white space Made widespread corrections to indentation and other white space issues to conform to the MOM6 standards, as documented in the white space section of the MOM6 style guide, at https://github.com/mom-ocean/MOM6/wiki/Code-style-guide. - MOM6 code uses a 2-point indent scheme. - There should be a space around all assignment (" = ") operators. - The names of optional arguments do not use whitespace, to differentiate them from assignments. - There is a space around the semicolons for stacked enddo statements. After this change, greping for '^ [a-zA-Z]', '^ [a-zA-Z]' or '^ [a-zA-Z]' does not find any lines in the MOM6 src or config_src/infra code. All answers and output are bitwise identical. --- .../infra/FMS1/MOM_diag_manager_infra.F90 | 60 +-- config_src/infra/FMS1/MOM_domain_infra.F90 | 16 +- .../infra/FMS2/MOM_diag_manager_infra.F90 | 60 +-- config_src/infra/FMS2/MOM_domain_infra.F90 | 16 +- src/ALE/MOM_ALE.F90 | 2 +- src/ALE/MOM_hybgen_remap.F90 | 2 +- src/ALE/MOM_remapping.F90 | 2 +- src/core/MOM.F90 | 2 +- src/core/MOM_barotropic.F90 | 2 +- src/core/MOM_dynamics_split_RK2.F90 | 4 +- src/core/MOM_dynamics_unsplit.F90 | 4 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 4 +- src/core/MOM_forcing_type.F90 | 14 +- src/core/MOM_open_boundary.F90 | 2 +- src/diagnostics/MOM_diagnostics.F90 | 2 +- src/equation_of_state/MOM_EOS_NEMO.F90 | 4 +- src/equation_of_state/MOM_TFreeze.F90 | 2 +- src/framework/MOM_diag_mediator.F90 | 2 +- src/framework/MOM_domains.F90 | 2 +- src/framework/MOM_horizontal_regridding.F90 | 2 +- src/ice_shelf/MOM_ice_shelf.F90 | 20 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 198 ++++----- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 342 +++++++------- .../MOM_state_initialization.F90 | 28 +- src/ocean_data_assim/MOM_oda_incupd.F90 | 420 +++++++++--------- .../vertical/MOM_ALE_sponge.F90 | 6 +- .../vertical/MOM_CVMix_conv.F90 | 2 +- .../vertical/MOM_CVMix_ddiff.F90 | 2 +- .../vertical/MOM_CVMix_shear.F90 | 4 +- .../vertical/MOM_bkgnd_mixing.F90 | 2 +- .../vertical/MOM_diabatic_driver.F90 | 2 +- src/tracer/MOM_OCMIP2_CFC.F90 | 10 +- src/tracer/MOM_generic_tracer.F90 | 192 ++++---- src/tracer/MOM_lateral_boundary_diffusion.F90 | 12 +- src/tracer/MOM_neutral_diffusion.F90 | 4 +- src/tracer/MOM_tracer_hor_diff.F90 | 2 +- src/user/BFB_surface_forcing.F90 | 2 +- src/user/dumbbell_surface_forcing.F90 | 2 +- 38 files changed, 726 insertions(+), 728 deletions(-) diff --git a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 index 18c80cf24c..69d82b7d85 100644 --- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 @@ -72,36 +72,36 @@ module MOM_diag_manager_infra !> Initialize a diagnostic axis integer function MOM_diag_axis_init(name, data, units, cart_name, long_name, MOM_domain, position, & & direction, edges, set_name, coarsen, null_axis) - character(len=*), intent(in) :: name !< The name of this axis - real, dimension(:), intent(in) :: data !< The array of coordinate values - character(len=*), intent(in) :: units !< The units for the axis data - character(len=*), intent(in) :: cart_name !< Cartesian axis ("X", "Y", "Z", "T", or "N" for none) - character(len=*), & - optional, intent(in) :: long_name !< The long name of this axis - type(MOM_domain_type), & - optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition - integer, optional, intent(in) :: position !< This indicates the relative position of this - !! axis. The default is CENTER, but EAST and NORTH - !! are common options. - integer, optional, intent(in) :: direction !< This indicates the direction along which this - !! axis increases: 1 for upward, -1 for downward, or - !! 0 for non-vertical axes (the default) - integer, optional, intent(in) :: edges !< The axis_id of the complementary axis that - !! describes the edges of this axis - character(len=*), & - optional, intent(in) :: set_name !< A name to use for this set of axes. - integer, optional, intent(in) :: coarsen !< An optional degree of coarsening for the grid, 1 - !! by default. - logical, optional, intent(in) :: null_axis !< If present and true, return the special null axis - !! id for use with scalars. - - integer :: coarsening ! The degree of grid coarsening - - if (present(null_axis)) then ; if (null_axis) then - ! Return the special null axis id for scalars - MOM_diag_axis_init = null_axis_id - return - endif ; endif + character(len=*), intent(in) :: name !< The name of this axis + real, dimension(:), intent(in) :: data !< The array of coordinate values + character(len=*), intent(in) :: units !< The units for the axis data + character(len=*), intent(in) :: cart_name !< Cartesian axis ("X", "Y", "Z", "T", or "N" for none) + character(len=*), & + optional, intent(in) :: long_name !< The long name of this axis + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition + integer, optional, intent(in) :: position !< This indicates the relative position of this + !! axis. The default is CENTER, but EAST and NORTH + !! are common options. + integer, optional, intent(in) :: direction !< This indicates the direction along which this + !! axis increases: 1 for upward, -1 for downward, or + !! 0 for non-vertical axes (the default) + integer, optional, intent(in) :: edges !< The axis_id of the complementary axis that + !! describes the edges of this axis + character(len=*), & + optional, intent(in) :: set_name !< A name to use for this set of axes. + integer, optional, intent(in) :: coarsen !< An optional degree of coarsening for the grid, 1 + !! by default. + logical, optional, intent(in) :: null_axis !< If present and true, return the special null axis + !! id for use with scalars. + + integer :: coarsening ! The degree of grid coarsening + + if (present(null_axis)) then ; if (null_axis) then + ! Return the special null axis id for scalars + MOM_diag_axis_init = null_axis_id + return + endif ; endif if (present(MOM_domain)) then coarsening = 1 ; if (present(coarsen)) coarsening = coarsen diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90 index f28c36b9d2..249453f42b 100644 --- a/config_src/infra/FMS1/MOM_domain_infra.F90 +++ b/config_src/infra/FMS1/MOM_domain_infra.F90 @@ -1906,14 +1906,14 @@ end subroutine get_simple_array_j_ind !> Invert the contents of a 1-d array subroutine invert(array) - integer, dimension(:), intent(inout) :: array !< The 1-d array to invert - integer :: i, ni, swap - ni = size(array) - do i=1,ni - swap = array(i) - array(i) = array(ni+1-i) - array(ni+1-i) = swap - enddo + integer, dimension(:), intent(inout) :: array !< The 1-d array to invert + integer :: i, ni, swap + ni = size(array) + do i=1,ni + swap = array(i) + array(i) = array(ni+1-i) + array(ni+1-i) = swap + enddo end subroutine invert !> Returns the global shape of h-point arrays diff --git a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 index 18c80cf24c..69d82b7d85 100644 --- a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 @@ -72,36 +72,36 @@ module MOM_diag_manager_infra !> Initialize a diagnostic axis integer function MOM_diag_axis_init(name, data, units, cart_name, long_name, MOM_domain, position, & & direction, edges, set_name, coarsen, null_axis) - character(len=*), intent(in) :: name !< The name of this axis - real, dimension(:), intent(in) :: data !< The array of coordinate values - character(len=*), intent(in) :: units !< The units for the axis data - character(len=*), intent(in) :: cart_name !< Cartesian axis ("X", "Y", "Z", "T", or "N" for none) - character(len=*), & - optional, intent(in) :: long_name !< The long name of this axis - type(MOM_domain_type), & - optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition - integer, optional, intent(in) :: position !< This indicates the relative position of this - !! axis. The default is CENTER, but EAST and NORTH - !! are common options. - integer, optional, intent(in) :: direction !< This indicates the direction along which this - !! axis increases: 1 for upward, -1 for downward, or - !! 0 for non-vertical axes (the default) - integer, optional, intent(in) :: edges !< The axis_id of the complementary axis that - !! describes the edges of this axis - character(len=*), & - optional, intent(in) :: set_name !< A name to use for this set of axes. - integer, optional, intent(in) :: coarsen !< An optional degree of coarsening for the grid, 1 - !! by default. - logical, optional, intent(in) :: null_axis !< If present and true, return the special null axis - !! id for use with scalars. - - integer :: coarsening ! The degree of grid coarsening - - if (present(null_axis)) then ; if (null_axis) then - ! Return the special null axis id for scalars - MOM_diag_axis_init = null_axis_id - return - endif ; endif + character(len=*), intent(in) :: name !< The name of this axis + real, dimension(:), intent(in) :: data !< The array of coordinate values + character(len=*), intent(in) :: units !< The units for the axis data + character(len=*), intent(in) :: cart_name !< Cartesian axis ("X", "Y", "Z", "T", or "N" for none) + character(len=*), & + optional, intent(in) :: long_name !< The long name of this axis + type(MOM_domain_type), & + optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition + integer, optional, intent(in) :: position !< This indicates the relative position of this + !! axis. The default is CENTER, but EAST and NORTH + !! are common options. + integer, optional, intent(in) :: direction !< This indicates the direction along which this + !! axis increases: 1 for upward, -1 for downward, or + !! 0 for non-vertical axes (the default) + integer, optional, intent(in) :: edges !< The axis_id of the complementary axis that + !! describes the edges of this axis + character(len=*), & + optional, intent(in) :: set_name !< A name to use for this set of axes. + integer, optional, intent(in) :: coarsen !< An optional degree of coarsening for the grid, 1 + !! by default. + logical, optional, intent(in) :: null_axis !< If present and true, return the special null axis + !! id for use with scalars. + + integer :: coarsening ! The degree of grid coarsening + + if (present(null_axis)) then ; if (null_axis) then + ! Return the special null axis id for scalars + MOM_diag_axis_init = null_axis_id + return + endif ; endif if (present(MOM_domain)) then coarsening = 1 ; if (present(coarsen)) coarsening = coarsen diff --git a/config_src/infra/FMS2/MOM_domain_infra.F90 b/config_src/infra/FMS2/MOM_domain_infra.F90 index a19f022cc7..d845d7317b 100644 --- a/config_src/infra/FMS2/MOM_domain_infra.F90 +++ b/config_src/infra/FMS2/MOM_domain_infra.F90 @@ -1905,14 +1905,14 @@ end subroutine get_simple_array_j_ind !> Invert the contents of a 1-d array subroutine invert(array) - integer, dimension(:), intent(inout) :: array !< The 1-d array to invert - integer :: i, ni, swap - ni = size(array) - do i=1,ni - swap = array(i) - array(i) = array(ni+1-i) - array(ni+1-i) = swap - enddo + integer, dimension(:), intent(inout) :: array !< The 1-d array to invert + integer :: i, ni, swap + ni = size(array) + do i=1,ni + swap = array(i) + array(i) = array(ni+1-i) + array(ni+1-i) = swap + enddo end subroutine invert !> Returns the global shape of h-point arrays diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 876c2c684b..e5666a5157 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -334,7 +334,7 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) 'm', conversion=GV%H_to_m, v_extensive=.true.) cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, & 'Layer thicknesses tendency due to ALE regridding and remapping', & - 'm s-1', conversion=GV%H_to_m*US%s_to_T, v_extensive = .true.) + 'm s-1', conversion=GV%H_to_m*US%s_to_T, v_extensive=.true.) end subroutine ALE_register_diags diff --git a/src/ALE/MOM_hybgen_remap.F90 b/src/ALE/MOM_hybgen_remap.F90 index 539c992283..59974afad1 100644 --- a/src/ALE/MOM_hybgen_remap.F90 +++ b/src/ALE/MOM_hybgen_remap.F90 @@ -192,7 +192,7 @@ subroutine hybgen_ppm_coefs(s, h_src, edges, nk, ns, thin, PCM_lay) + I_h0123(K)*( 2.*dp(k)*dp(k-1)*I_h12(K)*(s(k,i)-s(k-1,i)) * & ( h01_h112(K) - h23_h122(K) ) & + (dp(k)*as(k-1)*h23_h122(K) - dp(k-1)*as(k)*h01_h112(K)) ) - ar(k-1) = al(k) + ar(k-1) = al(k) enddo !k ar(nk-1) = s(nk,i) ! last layer PCM al(nk) = s(nk,i) ! last layer PCM diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index ebe0f81743..4146237e21 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -205,7 +205,7 @@ subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edg real, optional, intent(in) :: h_neglect_edge !< A negligibly small width !! for the purpose of edge value !! calculations in the same units as h0 [H] - logical, dimension(n0), optional, intent(in) :: PCM_cell !< If present, use PCM remapping for + logical, dimension(n0), optional, intent(in) :: PCM_cell !< If present, use PCM remapping for !! cells in the source grid where this is true. ! Local variables diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0703e45696..db08802cca 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1908,7 +1908,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "If False, T/S are registered for advection. "//& "This is intended only to be used in offline tracer mode "//& "and is by default false in that case.", & - do_not_log = .true., default=.true. ) + do_not_log=.true., default=.true.) if (present(offline_tracer_mode)) then ! Only read this parameter in enabled modes call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", CS%offline_tracer_mode, & "If true, barotropic and baroclinic dynamics, thermodynamics "//& diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 3b5c812ba1..d4e365f688 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -3053,7 +3053,7 @@ end subroutine apply_velocity_OBCs !! boundary conditions, as developed by Mehmet Ilicak. subroutine set_up_BT_OBC(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_BT_cont, & integral_BT_cont, dt_baroclinic, Datu, Datv, BTCL_u, BTCL_v) - type(ocean_OBC_type), target, intent(inout) :: OBC !< An associated pointer to an OBC type. + type(ocean_OBC_type), target, intent(inout) :: OBC !< An associated pointer to an OBC type. type(memory_size_type), intent(in) :: MS !< A type that describes the memory sizes of the !! argument arrays. real, dimension(SZIW_(MS),SZJW_(MS)), intent(in) :: eta !< The barotropic free surface height anomaly or diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index d177e63be8..9f13268090 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1330,10 +1330,10 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - y_cell_method='sum', v_extensive = .true.) + y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - x_cell_method='sum', v_extensive = .true.) + x_cell_method='sum', v_extensive=.true.) !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, & ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 085b4e02c3..e23a2c4ca1 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -694,10 +694,10 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & 'Effective U Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - y_cell_method='sum', v_extensive = .true.) + y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & 'Effective V Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - x_cell_method='sum', v_extensive = .true.) + x_cell_method='sum', v_extensive=.true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 91aaca508e..8c855d3aec 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -656,10 +656,10 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, & 'Effective U-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - y_cell_method='sum', v_extensive = .true.) + y_cell_method='sum', v_extensive=.true.) CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, & 'Effective V-Face Area', 'm^2', conversion = GV%H_to_m*US%L_to_m, & - x_cell_method='sum', v_extensive = .true.) + x_cell_method='sum', v_extensive=.true.) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 6de46f7738..9fe449e726 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -265,15 +265,15 @@ module MOM_forcing_type !! ice needs to be accumulated, and the rigidity explicitly !! reset to zero at the driver level when appropriate. real, pointer, dimension(:,:) :: & - ustk0 => NULL(), & !< Surface Stokes drift, zonal [m s-1] - vstk0 => NULL() !< Surface Stokes drift, meridional [m s-1] + ustk0 => NULL(), & !< Surface Stokes drift, zonal [m s-1] + vstk0 => NULL() !< Surface Stokes drift, meridional [m s-1] real, pointer, dimension(:) :: & - stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad m-1] + stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad m-1] real, pointer, dimension(:,:,:) :: & - ustkb => NULL(), & !< Stokes Drift spectrum, zonal [m s-1] + ustkb => NULL(), & !< Stokes Drift spectrum, zonal [m s-1] !! Horizontal - u points !! 3rd dimension - wavenumber - vstkb => NULL() !< Stokes Drift spectrum, meridional [m s-1] + vstkb => NULL() !< Stokes Drift spectrum, meridional [m s-1] !! Horizontal - v points !! 3rd dimension - wavenumber @@ -1503,14 +1503,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, cmor_field_name='ave_evs', cmor_standard_name='water_evaporation_flux_area_averaged', & cmor_long_name='Evaporation Where Ice Free Ocean over Sea Area Averaged') - handles%id_lprec_ga = register_scalar_field('ocean_model', 'lprec_ga', Time, diag,& + handles%id_lprec_ga = register_scalar_field('ocean_model', 'lprec_ga', Time, diag,& long_name='Area integrated liquid precip into ocean', & units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & standard_name='rainfall_flux_area_averaged', & cmor_field_name='ave_pr', cmor_standard_name='rainfall_flux_area_averaged', & cmor_long_name='Rainfall Flux where Ice Free Ocean over Sea Area Averaged') - handles%id_fprec_ga = register_scalar_field('ocean_model', 'fprec_ga', Time, diag, & + handles%id_fprec_ga = register_scalar_field('ocean_model', 'fprec_ga', Time, diag, & long_name='Area integrated frozen precip into ocean', & units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s, & standard_name='snowfall_flux_area_averaged', & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 39ceccaf49..971fd7e8ee 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1630,7 +1630,7 @@ end subroutine parse_segment_data_str !> Parse all the OBC_SEGMENT_%%%_DATA strings again !! to see which need tracer reservoirs (all pes need to know). - subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) +subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) type(ocean_OBC_type), target, intent(inout) :: OBC !< Open boundary control structure type(param_file_type), intent(in) :: PF !< Parameter file handle logical, intent(in) :: use_temperature !< If true, T and S are used diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 22c0ddf8de..9c22ab5eb5 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -452,7 +452,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! area mean SSS if (CS%id_sosga > 0) then do j=js,je ; do i=is,ie - surface_field(i,j) = tv%S(i,j,1) + surface_field(i,j) = tv%S(i,j,1) enddo ; enddo sosga = global_area_mean(surface_field, G) call post_data(CS%id_sosga, sosga, CS%diag) diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index 68488881bb..b8c2ec2601 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -259,7 +259,7 @@ subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts, rho_re rho(j) = ( zn + zr0 ) ! density endif - enddo + enddo end subroutine calculate_density_array_nemo !> For a given thermodynamic state, calculate the derivatives of density with conservative @@ -391,7 +391,7 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) - enddo + enddo end subroutine calculate_compress_nemo end module MOM_EOS_NEMO diff --git a/src/equation_of_state/MOM_TFreeze.F90 b/src/equation_of_state/MOM_TFreeze.F90 index 50233cae60..f0b22c8f4e 100644 --- a/src/equation_of_state/MOM_TFreeze.F90 +++ b/src/equation_of_state/MOM_TFreeze.F90 @@ -165,7 +165,7 @@ subroutine calculate_TFreeze_teos10_array(S, pres, T_Fr, start, npts) if (S(j) < -1.0e-10) cycle !Can we assume safely that this is a missing value? T_Fr(j) = gsw_ct_freezing_exact(zs,zp,saturation_fraction) - enddo + enddo end subroutine calculate_TFreeze_teos10_array diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 7fb2580906..be8ccd774b 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -2234,7 +2234,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time ! Register the native diagnostic if (associated(axes_d2)) then - active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, & + active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index a97ae52e35..47ac43df06 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -182,7 +182,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & !$ "The number of OpenMP threads that MOM6 will use.", & !$ default = 1, layoutParam=.true.) !$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, & - !$ "If True, use hyper-threading.", default = .false., layoutParam=.true.) + !$ "If True, use hyper-threading.", default=.false., layoutParam=.true.) !$ call set_MOM_thread_affinity(ocean_nthreads, ocean_omp_hyper_thread) !$ endif # endif diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 3bc16a6aff..1889c59469 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -249,7 +249,7 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, call MOM_error(FATAL,"MOM_initialize: "// & "fill is true and good is false after fill_miss, how did this happen? ") endif - enddo ; enddo + enddo ; enddo end subroutine fill_miss_2d diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 5db9198a22..5b63851b8e 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -327,8 +327,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) ISS => CS%ISS if (CS%data_override_shelf_fluxes .and. CS%active_shelf_dynamics) then - call data_override(G%Domain, 'shelf_sfc_mass_flux', fluxes_in%shelf_sfc_mass_flux, CS%Time, & - scale=US%kg_m2s_to_RZ_T) + call data_override(G%Domain, 'shelf_sfc_mass_flux', fluxes_in%shelf_sfc_mass_flux, CS%Time, & + scale=US%kg_m2s_to_RZ_T) endif if (CS%rotate_index) then @@ -741,13 +741,13 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) scale=US%RZ_to_kg_m2) endif - call change_thickness_using_precip(CS, ISS, G, US, fluxes,US%s_to_T*time_step, Time) + call change_thickness_using_precip(CS, ISS, G, US, fluxes, US%s_to_T*time_step, Time) - if (CS%debug) then - call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, scale=US%Z_to_m) - call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using surf acc", G%HI, haloshift=0, & - scale=US%RZ_to_kg_m2) - endif + if (CS%debug) then + call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, scale=US%Z_to_m) + call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using surf acc", G%HI, haloshift=0, & + scale=US%RZ_to_kg_m2) + endif endif @@ -2136,8 +2136,8 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_ endif if (present(data_override_shelf_fluxes)) then - data_override_shelf_fluxes=.false. - if (CS%active_shelf_dynamics) data_override_shelf_fluxes = CS%data_override_shelf_fluxes + data_override_shelf_fluxes=.false. + if (CS%active_shelf_dynamics) data_override_shelf_fluxes = CS%data_override_shelf_fluxes endif end subroutine ice_shelf_query diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index e1ecd45347..979517d227 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -540,31 +540,31 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif ! initialize basal friction coefficients - if (new_sim) then - call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) - call pass_var(CS%C_basal_friction, G%domain) - - ! initialize ice-stiffness AGlen - call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) - call pass_var(CS%AGlen_visc, G%domain) - - !initialize boundary conditions - call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & - CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, & - CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file ) - call pass_var(ISS%hmask, G%domain) - call pass_var(CS%h_bdry_val, G%domain) - call pass_var(CS%thickness_bdry_val, G%domain) - call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) - call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) - - !initialize ice flow characteristic (velocities, bed elevation under the grounded part, etc) from file - call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, & - G, US, param_file) - call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) - call pass_var(CS%bed_elev, G%domain,CENTER) - call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) - endif + if (new_sim) then + call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) + call pass_var(CS%C_basal_friction, G%domain) + + ! initialize ice-stiffness AGlen + call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) + call pass_var(CS%AGlen_visc, G%domain) + + !initialize boundary conditions + call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & + CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, & + CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file ) + call pass_var(ISS%hmask, G%domain) + call pass_var(CS%h_bdry_val, G%domain) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) + call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) + + !initialize ice flow characteristic (velocities, bed elevation under the grounded part, etc) from file + call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf, CS%ground_frac, & + G, US, param_file) + call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) + call pass_var(CS%bed_elev, G%domain,CENTER) + call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) + endif ! Register diagnostics. CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesB1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) @@ -628,7 +628,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) enddo enddo - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) end subroutine initialize_diagnostic_fields !> This function returns the global maximum advective timestep that can be taken based on the current @@ -691,58 +691,58 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding ! - call ice_shelf_advect(CS, ISS, G, time_step, Time) - CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step - if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. - - if (coupled_GL) then - call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel) - elseif (update_ice_vel) then - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - endif - - - if (update_ice_vel) then - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) - endif - -! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) - - if (update_ice_vel) then - call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) - if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) - if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) - if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) -! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) - if (CS%id_taudx_shelf > 0) then - taud_x(:,:) = CS%taudx_shelf(:,:)*G%IareaT(:,:) - call post_data(CS%id_taudx_shelf,taud_x , CS%diag) - endif - if (CS%id_taudy_shelf > 0) then - taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaT(:,:) - call post_data(CS%id_taudy_shelf,taud_y , CS%diag) - endif - if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) - if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) - if (CS%id_visc_shelf > 0) then - ice_visc(:,:)=CS%ice_visc(:,:)*G%IareaT(:,:) - call post_data(CS%id_visc_shelf, ice_visc,CS%diag) - endif - if (CS%id_taub > 0) then - basal_tr(:,:) = CS%basal_traction(:,:)*G%IareaT(:,:) - call post_data(CS%id_taub, basal_tr,CS%diag) - endif + call ice_shelf_advect(CS, ISS, G, time_step, Time) + CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step + if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. + + if (coupled_GL) then + call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel) + elseif (update_ice_vel) then + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + endif + + + if (update_ice_vel) then + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) + endif + +! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) + + if (update_ice_vel) then + call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) + if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) + if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) + if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) +! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf, CS%t_shelf, CS%diag) + if (CS%id_taudx_shelf > 0) then + taud_x(:,:) = CS%taudx_shelf(:,:)*G%IareaT(:,:) + call post_data(CS%id_taudx_shelf, taud_x, CS%diag) + endif + if (CS%id_taudy_shelf > 0) then + taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaT(:,:) + call post_data(CS%id_taudy_shelf, taud_y, CS%diag) + endif + if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac, CS%diag) + if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) + if (CS%id_visc_shelf > 0) then + ice_visc(:,:) = CS%ice_visc(:,:)*G%IareaT(:,:) + call post_data(CS%id_visc_shelf, ice_visc, CS%diag) + endif + if (CS%id_taub > 0) then + basal_tr(:,:) = CS%basal_traction(:,:)*G%IareaT(:,:) + call post_data(CS%id_taub, basal_tr, CS%diag) + endif !! - if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) - if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) - if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag) - if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag) -! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask,CS%tmask,CS%diag) + if (CS%id_u_mask > 0) call post_data(CS%id_u_mask, CS%umask, CS%diag) + if (CS%id_v_mask > 0) call post_data(CS%id_v_mask, CS%vmask, CS%diag) + if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask, CS%u_face_mask_bdry, CS%diag) + if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask, CS%v_face_mask_bdry, CS%diag) +! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask, CS%tmask, CS%diag) - call disable_averaging(CS%diag) + call disable_averaging(CS%diag) - CS%elapsed_velocity_time = 0.0 - endif + CS%elapsed_velocity_time = 0.0 + endif end subroutine update_ice_shelf @@ -836,7 +836,7 @@ end subroutine ice_shelf_advect !>This subroutine computes u- and v-velocities of the ice shelf iterating on non-linear ice viscosity !subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) - subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, iters, Time) +subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, iters, Time) type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state @@ -948,7 +948,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i ! This makes sure basal stress is only applied when it is supposed to be do j=G%jsd,G%jed ; do i=G%isd,G%ied ! CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) + CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) enddo ; enddo call apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, CS%ice_visc, & @@ -1001,8 +1001,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i ! makes sure basal stress is only applied when it is supposed to be do j=G%jsd,G%jed ; do i=G%isd,G%ied -! CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) +! CS%basal_traction(i,j) = CS%basal_traction(i,j) * CS%ground_frac(i,j) + CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) enddo ; enddo u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0 @@ -1371,7 +1371,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H enddo call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) - if (conv_flag == 0) then + if (conv_flag == 0) then iters = CS%cg_max_iterations endif @@ -1859,8 +1859,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo enddo - call pass_var(S, G%domain) - allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) + call pass_var(S, G%domain) + allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) do j=jscq,jecq ; do i=iscq,iecq call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j)) enddo ; enddo @@ -2011,7 +2011,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo enddo - deallocate(Phi) + deallocate(Phi) end subroutine calc_shelf_driving_stress subroutine init_boundary_values(CS, G, time, hmask, input_flux, input_thick, new_sim) @@ -2616,11 +2616,10 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) n_g = CS%n_glen; eps_min = CS%eps_glen_min CS%ice_visc(:,:)=1e22 ! 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 + do j=jsc,jec ; do i=isc,iec - if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then - Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%AGlen_visc(i,j))**(-1./CS%n_glen) + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then + Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%AGlen_visc(i,j))**(-1./CS%n_glen) do iq=1,2 ; do jq=1,2 @@ -2643,13 +2642,12 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j)) + & (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 -! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging - CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + enddo ; enddo +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging + 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)) - endif - enddo - enddo + endif + enddo ; enddo deallocate(Phi) end subroutine calc_shelf_visc @@ -2863,16 +2861,16 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3)) do qpoint=1,4 - if (J>1) then + if (J>1) then a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) - else - a= G%dxCv(i,J) !* yquad(qpoint) ! d(x)/d(x*) - endif - if (I>1) then + else + a = G%dxCv(i,J) !* yquad(qpoint) ! d(x)/d(x*) + endif + if (I>1) then d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) - else + else d = G%dyCu(I,j) !* xquad(qpoint) - endif + endif ! a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) ! d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index d3202a9f2a..7bd3b45b2a 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -277,120 +277,120 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b thickness_bdry_val, hmask, h_shelf, G,& US, PF ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces - - real, dimension(SZIB_(G),SZJ_(G)), & - intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through - !! C-grid u faces [L Z T-1 ~> m2 s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces - - real, dimension(SZI_(G),SZJB_(G)), & - intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through - !! C-grid v faces [L Z T-1 ~> m2 s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] - !! boundary vertices [L T-1 ~> m s-1]. - - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_shelf !< Ice-shelf thickness [Z ~> m] - type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors - type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters - - character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. - integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed - real :: input_thick ! The input ice shelf thickness [Z ~> m] - real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] - real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises - - - call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) - - call get_param(PF, mdl, "LENLON", lenlon, fail_if_missing=.true.) - - call get_param(PF, mdl, "WESTLON", westlon, fail_if_missing=.true.) - - call get_param(PF, mdl, "SOUTHLAT", southlat, fail_if_missing=.true.) - - call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & - "inflow ice velocity at upstream boundary", & - units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) !### This conversion factor is wrong? - call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & - "flux thickness at upstream boundary", & - units="m", default=1000., scale=US%m_to_Z) - call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & - "maximum position of no-flow condition in along-flow direction", & - units="km", default=0.) - - call MOM_mesg(mdl//": setting boundary") - - isd = G%isd ; ied = G%ied - jsd = G%jsd ; jed = G%jed - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec - gjsd = G%Domain%njglobal ; gisd = G%Domain%niglobal - gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo - giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc - - !---------b.c.s based on geopositions ----------------- - do j=jsc,jec+1 - do i=isc-1,iec+1 - ! upstream boundary - set either dirichlet or flux condition - - if (G%geoLonBu(i,j) == westlon) then - hmask(i+1,j) = 3.0 - h_bdry_val(i+1,j) = h_shelf(i+1,j) - thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) - u_face_mask_bdry(i+1,j) = 3.0 - u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution - endif - - - ! side boundaries: no flow - if (G%geoLatBu(i,j-1) == southlat) then !bot boundary - if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then - v_face_mask_bdry(i,j+1) = 0. - u_face_mask_bdry(i,j) = 3. - u_bdry_val(i,j) = 0. - v_bdry_val(i,j) = 0. - else - v_face_mask_bdry(i,j+1) = 1. - u_face_mask_bdry(i,j) = 3. - u_bdry_val(i,j) = 0. - v_bdry_val(i,j) = 0. - endif - elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary - if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then - v_face_mask_bdry(i,j-1) = 0. - u_face_mask_bdry(i,j-1) = 3. - else - v_face_mask_bdry(i,j-1) = 3. - u_face_mask_bdry(i,j-1) = 3. - endif - endif - - ! downstream boundary - CFBC - if (G%geoLonBu(i,j) == westlon+lenlon) then - u_face_mask_bdry(i-1,j) = 2.0 - endif - - enddo - enddo + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces + + real, dimension(SZIB_(G),SZJ_(G)), & + intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through + !! C-grid u faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces + + real, dimension(SZI_(G),SZJB_(G)), & + intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through + !! C-grid v faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + !! boundary vertices [L T-1 ~> m s-1]. + + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_shelf !< Ice-shelf thickness [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. + integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed + real :: input_thick ! The input ice shelf thickness [Z ~> m] + real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] + real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises + + + call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "LENLON", lenlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "WESTLON", westlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "SOUTHLAT", southlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & + "inflow ice velocity at upstream boundary", & + units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) !### This conversion factor is wrong? + call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & + "flux thickness at upstream boundary", & + units="m", default=1000., scale=US%m_to_Z) + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & + "maximum position of no-flow condition in along-flow direction", & + units="km", default=0.) + + call MOM_mesg(mdl//": setting boundary") + + isd = G%isd ; ied = G%ied + jsd = G%jsd ; jed = G%jed + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + gjsd = G%Domain%njglobal ; gisd = G%Domain%niglobal + gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo + giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc + + !---------b.c.s based on geopositions ----------------- + do j=jsc,jec+1 + do i=isc-1,iec+1 + ! upstream boundary - set either dirichlet or flux condition + + if (G%geoLonBu(i,j) == westlon) then + hmask(i+1,j) = 3.0 + h_bdry_val(i+1,j) = h_shelf(i+1,j) + thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) + u_face_mask_bdry(i+1,j) = 3.0 + u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution + endif + + + ! side boundaries: no flow + if (G%geoLatBu(i,j-1) == southlat) then !bot boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j+1) = 0. + u_face_mask_bdry(i,j) = 3. + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. + else + v_face_mask_bdry(i,j+1) = 1. + u_face_mask_bdry(i,j) = 3. + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. + endif + elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j-1) = 0. + u_face_mask_bdry(i,j-1) = 3. + else + v_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j-1) = 3. + endif + endif + + ! downstream boundary - CFBC + if (G%geoLonBu(i,j) == westlon+lenlon) then + u_face_mask_bdry(i-1,j) = 2.0 + endif + + enddo + enddo end subroutine initialize_ice_shelf_boundary_channel @@ -400,10 +400,10 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: bed_elev !< The bed elevation [Z ~> m]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: float_cond !< An array indicating where the ice @@ -471,33 +471,33 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask u_bdry_val, v_bdry_val, umask, vmask, h_bdry_val, thickness_bdry_val, & hmask, h_shelf, G, US, PF ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_face_mask_bdry !< A boundary-type mask at B-grid u faces [nondim] - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_face_mask_bdry !< A boundary-type mask at B-grid v faces [nondim] - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. - real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. - real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: umask !< A mask for ice shelf velocity [nondim] - real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: vmask !< A mask for ice shelf velocity [nondim] - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] - !! boundary vertices [L T-1 ~> m s-1]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf [nondim] - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: h_shelf !< Ice-shelf thickness [Z ~> m] - type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors - type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_face_mask_bdry !< A boundary-type mask at B-grid u faces [nondim] + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_face_mask_bdry !< A boundary-type mask at B-grid v faces [nondim] + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(inout) :: umask !< A mask for ice shelf velocity [nondim] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(inout) :: vmask !< A mask for ice shelf velocity [nondim] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf [nondim] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: h_shelf !< Ice-shelf thickness [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters character(len=200) :: filename, bc_file, inputdir, icethick_file ! Strings for file/path character(len=200) :: ufcmskbdry_varname, vfcmskbdry_varname, & @@ -562,10 +562,10 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask do j=jsc,jec do i=isc,iec - if (hmask(i,j) == 3.) then - thickness_bdry_val(i,j) = h_shelf(i,j) - h_bdry_val(i,j) = h_shelf(i,j) - endif + if (hmask(i,j) == 3.) then + thickness_bdry_val(i,j) = h_shelf(i,j) + h_bdry_val(i,j) = h_shelf(i,j) + endif enddo enddo @@ -574,8 +574,8 @@ end subroutine initialize_ice_shelf_boundary_from_file !> Initialize ice basal friction subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: C_basal_friction !< Ice-stream basal friction + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: C_basal_friction !< Ice-stream basal friction type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -595,19 +595,19 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) call get_param(PF, mdl, "BASAL_FRICTION_COEFF", C_friction, & "Coefficient in sliding law.", units="Pa (m s-1)^(n_basal_fric)", default=5.e10) - C_basal_friction(:,:) = C_friction + C_basal_friction(:,:) = C_friction elseif (trim(config)=="FILE") then - call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading friction coefficients") - call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) + call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading friction coefficients") + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) - call get_param(PF, mdl, "BASAL_FRICTION_FILE", C_friction_file, & - "The file from which basal friction coefficients are read.", & - default="ice_basal_friction.nc") - filename = trim(inputdir)//trim(C_friction_file) - call log_param(PF, mdl, "INPUTDIR/BASAL_FRICTION_FILE", filename) + call get_param(PF, mdl, "BASAL_FRICTION_FILE", C_friction_file, & + "The file from which basal friction coefficients are read.", & + default="ice_basal_friction.nc") + filename = trim(inputdir)//trim(C_friction_file) + call log_param(PF, mdl, "INPUTDIR/BASAL_FRICTION_FILE", filename) - call get_param(PF, mdl, "BASAL_FRICTION_VARNAME", varname, & + call get_param(PF, mdl, "BASAL_FRICTION_VARNAME", varname, & "The variable to use in basal traction.", & default="tau_b_beta") @@ -623,8 +623,8 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) !> Initialize ice basal friction subroutine initialize_ice_AGlen(AGlen, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: AGlen !< The ice-stiffness parameter A_Glen + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: AGlen !< The ice-stiffness parameter A_Glen type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -644,19 +644,19 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF) call get_param(PF, mdl, "A_GLEN", A_Glen, & "Ice-stiffness parameter.", units="Pa-3 s-1", default=2.261e-25) - AGlen(:,:) = A_Glen + AGlen(:,:) = A_Glen elseif (trim(config)=="FILE") then - call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading ice-stiffness parameter") - call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) + call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading ice-stiffness parameter") + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) - call get_param(PF, mdl, "ICE_STIFFNESS_FILE", AGlen_file, & + call get_param(PF, mdl, "ICE_STIFFNESS_FILE", AGlen_file, & "The file from which the ice-stiffness is read.", & default="ice_AGlen.nc") - filename = trim(inputdir)//trim(AGlen_file) - call log_param(PF, mdl, "INPUTDIR/ICE_STIFFNESS_FILE", filename) - call get_param(PF, mdl, "A_GLEN_VARNAME", varname, & + filename = trim(inputdir)//trim(AGlen_file) + call log_param(PF, mdl, "INPUTDIR/ICE_STIFFNESS_FILE", filename) + call get_param(PF, mdl, "A_GLEN_VARNAME", varname, & "The variable to use as ice-stiffness.", & default="A_GLEN") diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index ee3052f166..0aeba26e17 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1947,7 +1947,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t "The name of the inverse damping rate variable in "//& "SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var) endif - call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log=.true.) !### NEW_SPONGES should be obsoleted properly, rather than merely deprecated, at which ! point only the else branch of the new_sponge_param block would be retained. @@ -2000,18 +2000,18 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t call MOM_read_vector(filename, Idamp_u_var,Idamp_v_var,Idamp_u(:,:),Idamp_v(:,:), G%Domain, scale=US%T_to_s) else - ! call MOM_error(FATAL, "Must provide SPONGE_IDAMP_U_var and SPONGE_IDAMP_V_var") - call pass_var(Idamp,G%Domain) - do j=G%jsc,G%jec - do i=G%iscB,G%iecB - Idamp_u(I,j) = 0.5*(Idamp(i,j)+Idamp(i+1,j)) - enddo - enddo - do j=G%jscB,G%jecB - do i=G%isc,G%iec - Idamp_v(i,J) = 0.5*(Idamp(i,j)+Idamp(i,j+1)) - enddo - enddo + ! call MOM_error(FATAL, "Must provide SPONGE_IDAMP_U_var and SPONGE_IDAMP_V_var") + call pass_var(Idamp,G%Domain) + do j=G%jsc,G%jec + do i=G%iscB,G%iecB + Idamp_u(I,j) = 0.5*(Idamp(i,j)+Idamp(i+1,j)) + enddo + enddo + do j=G%jscB,G%jecB + do i=G%isc,G%iec + Idamp_v(i,J) = 0.5*(Idamp(i,j)+Idamp(i,j+1)) + enddo + enddo endif endif @@ -2255,7 +2255,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p "The name of the meridional vel. inc. variable in "//& "ODA_INCUPD_FILE.", default="v_inc") -! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) +! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log=.true.) ! Read in incremental update for tracers filename = trim(inputdir)//trim(inc_file) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 5f38fa15d0..80c313e488 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -201,13 +201,13 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re ! get number of timestep for full update if (nhours_incupd == 0) then - CS%nstep_incupd = 1 !! direct insertion + CS%nstep_incupd = 1 !! direct insertion else - CS%nstep_incupd = floor( nhours_incupd * 3600. / dt_therm + 0.001 ) - 1 + CS%nstep_incupd = floor( nhours_incupd * 3600. / dt_therm + 0.001 ) - 1 endif write(mesg,'(i12)') CS%nstep_incupd if (is_root_pe()) & - call MOM_error(NOTE,"initialize_oda_incupd: Number of Timestep of inc. update:"//& + call MOM_error(NOTE,"initialize_oda_incupd: Number of Timestep of inc. update:"//& trim(mesg)) ! number of inc. update already done, CS%ncount, either from restart or set to 0.0 @@ -219,14 +219,14 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re endif write(mesg,'(f4.1)') CS%ncount if (is_root_pe()) & - call MOM_error(NOTE,"initialize_oda_incupd: Inc. update already done:"//& + call MOM_error(NOTE,"initialize_oda_incupd: Inc. update already done:"//& trim(mesg)) ! get the vertical grid (h_obs) of the increments CS%nz_data = nz_data allocate(CS%Ref_h%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data), source=0.0) do j=G%jsc,G%jec; do i=G%isc,G%iec ; do k=1,CS%nz_data - CS%Ref_h%p(i,j,k) = data_h(i,j,k) + CS%Ref_h%p(i,j,k) = data_h(i,j,k) enddo; enddo ; enddo !### Doing a halo update here on CS%Ref_h%p would avoid needing halo updates each timestep. @@ -263,7 +263,7 @@ subroutine set_up_oda_incupd_field(sp_val, G, GV, CS) CS%Inc(CS%fldno)%nz_data = CS%nz_data allocate(CS%Inc(CS%fldno)%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data), source=0.0) do k=1,CS%nz_data ; do j=G%jsc,G%jec ; do i=G%isc,G%iec - CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) + CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) enddo ; enddo ; enddo end subroutine set_up_oda_incupd_field @@ -374,122 +374,122 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) ! remap t,s (on h_init) to h_obs to get increment tmp_val1(:) = 0.0 do j=js,je ; do i=is,ie - if (G%mask2dT(i,j) == 1) then + if (G%mask2dT(i,j) == 1) then + ! account for the different SSH + sum_h1 = 0.0 + sum_h2 = 0.0 + do k=1,nz + sum_h1 = sum_h1+h(i,j,k) + enddo + + do k=1,nz_data + sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data + tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) + enddo + ! get temperature + do k=1,nz + tmp_val1(k) = tv%T(i,j,k) + enddo + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & + nz_data, tmp_h(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc(1)%p(i,j,k) = CS%Inc(1)%p(i,j,k) - tmp_val2(k) + enddo + + ! get salinity + do k=1,nz + tmp_val1(k) = tv%S(i,j,k) + enddo + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & + nz_data, tmp_h(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc(2)%p(i,j,k) = CS%Inc(2)%p(i,j,k) - tmp_val2(k) + enddo + endif + enddo ; enddo + + ! remap u to h_obs to get increment + if (CS%uv_inc) then + call pass_var(h, G%Domain) + + hu(:) = 0.0 + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) == 1) then + ! get u-velocity + do k=1,nz + tmp_val1(k) = u(i,j,k) + ! get the h and h_obs at u points + hu(k) = 0.5*( h(i,j,k)+ h(i+1,j,k)) + enddo + do k=1,nz_data + hu_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i+1,j,k)) + enddo ! account for the different SSH sum_h1 = 0.0 - sum_h2 = 0.0 do k=1,nz - sum_h1 = sum_h1+h(i,j,k) + sum_h1 = sum_h1+hu(k) enddo - + sum_h2 = 0.0 do k=1,nz_data - sum_h2 = sum_h2+h_obs(i,j,k) + sum_h2 = sum_h2+hu_obs(k) enddo do k=1,nz_data - tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) + hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) enddo - ! get temperature - do k=1,nz - tmp_val1(k) = tv%T(i,j,k) - enddo - ! remap tracer on h_obs - call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & - nz_data, tmp_h(1:nz_data), tmp_val2, & + ! remap model u on hu_obs + call remapping_core_h(CS%remap_cs, nz, hu(1:nz), tmp_val1, & + nz_data, hu_obs(1:nz_data), tmp_val2, & h_neglect, h_neglect_edge) ! get increment from full field on h_obs do k=1,nz_data - CS%Inc(1)%p(i,j,k) = CS%Inc(1)%p(i,j,k) - tmp_val2(k) + CS%Inc_u%p(i,j,k) = CS%Inc_u%p(i,j,k) - tmp_val2(k) enddo - - ! get salinity + endif + enddo ; enddo + + ! remap v to h_obs to get increment + hv(:) = 0.0; + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) == 1) then + ! get v-velocity do k=1,nz - tmp_val1(k) = tv%S(i,j,k) + tmp_val1(k) = v(i,j,k) + ! get the h and h_obs at v points + hv(k) = 0.5*(h(i,j,k)+h(i,j+1,k)) enddo - ! remap tracer on h_obs - call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & - nz_data, tmp_h(1:nz_data), tmp_val2, & + do k=1,nz_data + hv_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i,j+1,k)) + enddo + ! account for the different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1+hv(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2+hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) + enddo + ! remap model v on hv_obs + call remapping_core_h(CS%remap_cs, nz, hv(1:nz), tmp_val1, & + nz_data, hv_obs(1:nz_data), tmp_val2, & h_neglect, h_neglect_edge) ! get increment from full field on h_obs do k=1,nz_data - CS%Inc(2)%p(i,j,k) = CS%Inc(2)%p(i,j,k) - tmp_val2(k) + CS%Inc_v%p(i,j,k) = CS%Inc_v%p(i,j,k) - tmp_val2(k) enddo - endif - enddo; enddo - - ! remap u to h_obs to get increment - if (CS%uv_inc) then - call pass_var(h, G%Domain) - - hu(:) = 0.0 - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(i,j) == 1) then - ! get u-velocity - do k=1,nz - tmp_val1(k) = u(i,j,k) - ! get the h and h_obs at u points - hu(k) = 0.5*( h(i,j,k)+ h(i+1,j,k)) - enddo - do k=1,nz_data - hu_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i+1,j,k)) - enddo - ! account for the different SSH - sum_h1 = 0.0 - do k=1,nz - sum_h1 = sum_h1+hu(k) - enddo - sum_h2 = 0.0 - do k=1,nz_data - sum_h2 = sum_h2+hu_obs(k) - enddo - do k=1,nz_data - hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) - enddo - ! remap model u on hu_obs - call remapping_core_h(CS%remap_cs, nz, hu(1:nz), tmp_val1, & - nz_data, hu_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - do k=1,nz_data - CS%Inc_u%p(i,j,k) = CS%Inc_u%p(i,j,k) - tmp_val2(k) - enddo - endif - enddo; enddo - - ! remap v to h_obs to get increment - hv(:) = 0.0; - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) == 1) then - ! get v-velocity - do k=1,nz - tmp_val1(k) = v(i,j,k) - ! get the h and h_obs at v points - hv(k) = 0.5*(h(i,j,k)+h(i,j+1,k)) - enddo - do k=1,nz_data - hv_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i,j+1,k)) - enddo - ! account for the different SSH - sum_h1 = 0.0 - do k=1,nz - sum_h1 = sum_h1+hv(k) - enddo - sum_h2 = 0.0 - do k=1,nz_data - sum_h2 = sum_h2+hv_obs(k) - enddo - do k=1,nz_data - hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) - enddo - ! remap model v on hv_obs - call remapping_core_h(CS%remap_cs, nz, hv(1:nz), tmp_val1, & - nz_data, hv_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - do k=1,nz_data - CS%Inc_v%p(i,j,k) = CS%Inc_v%p(i,j,k) - tmp_val2(k) - enddo - endif - enddo; enddo + endif + enddo ; enddo endif ! uv_inc call pass_var(CS%Inc(1)%p, G%Domain) @@ -603,118 +603,118 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) tmp_h(k) = ( sum_h1 / sum_h2 ) * h_obs(i,j,k) enddo if (G%mask2dT(i,j) == 1) then - ! get temperature increment - do k=1,nz_data - tmp_val2(k) = CS%Inc(1)%p(i,j,k) - enddo - ! remap increment profile on model h - call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, & - nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) - do k=1,nz - ! add increment to tracer on model h - tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k) - tmp_t(i,j,k) = tmp_val1(k) ! store T increment for diagnostics - enddo - - ! get salinity increment - do k=1,nz_data - tmp_val2(k) = CS%Inc(2)%p(i,j,k) - enddo - ! remap increment profile on model h - call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data),tmp_val2,& - nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) - ! add increment to tracer on model h - do k=1,nz - tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k) - tmp_s(i,j,k) = tmp_val1(k) ! store S increment for diagnostics - ! bound salinity values ! check if it is correct to do that or if it hides - ! other problems ... - tv%S(i,j,k) = max(0.0 , tv%S(i,j,k)) - enddo + ! get temperature increment + do k=1,nz_data + tmp_val2(k) = CS%Inc(1)%p(i,j,k) + enddo + ! remap increment profile on model h + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, & + nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) + do k=1,nz + ! add increment to tracer on model h + tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k) + tmp_t(i,j,k) = tmp_val1(k) ! store T increment for diagnostics + enddo + + ! get salinity increment + do k=1,nz_data + tmp_val2(k) = CS%Inc(2)%p(i,j,k) + enddo + ! remap increment profile on model h + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data),tmp_val2,& + nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) + ! add increment to tracer on model h + do k=1,nz + tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k) + tmp_s(i,j,k) = tmp_val1(k) ! store S increment for diagnostics + ! bound salinity values ! check if it is correct to do that or if it hides + ! other problems ... + tv%S(i,j,k) = max(0.0 , tv%S(i,j,k)) + enddo endif - enddo; enddo + enddo ; enddo ! add u and v increments if (CS%uv_inc) then - call pass_var(h,G%Domain) ! to ensure reproducibility - - ! add increments to u - hu(:) = 0.0 - tmp_u(:,:,:) = 0.0 ! diagnostics - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(i,j) == 1) then - do k=1,nz_data - ! get u increment - tmp_val2(k) = CS%Inc_u%p(i,j,k) - ! get the h and h_obs at u points - hu_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i+1,j,k) ) - enddo - do k=1,nz - hu(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k) ) - enddo - ! account for different SSH - sum_h1 = 0.0 - do k=1,nz - sum_h1 = sum_h1 + hu(k) - enddo - sum_h2 = 0.0 - do k=1,nz_data - sum_h2 = sum_h2 + hu_obs(k) - enddo - do k=1,nz_data - hu_obs(k)=( sum_h1 / sum_h2 ) * hu_obs(k) - enddo - ! remap increment profile on hu - call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & - nz, hu(1:nz), tmp_val1, h_neglect, h_neglect_edge) - ! add increment to u-velocity on hu - do k=1,nz - u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k) - ! store increment for diagnostics - tmp_u(i,j,k) = tmp_val1(k) - enddo - endif - enddo; enddo - - ! add increments to v - hv(:) = 0.0 - tmp_v(:,:,:) = 0.0 ! diagnostics - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) == 1) then - ! get v increment - do k=1,nz_data - tmp_val2(k) = CS%Inc_v%p(i,j,k) - ! get the h and h_obs at v points - hv_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i,j+1,k) ) - enddo - do k=1,nz - hv(k) = 0.5 * (h(i,j,k) + h(i,j+1,k) ) - enddo - ! account for different SSH - sum_h1 = 0.0 - do k=1,nz - sum_h1 = sum_h1 + hv(k) - enddo - sum_h2 = 0.0 - do k=1,nz_data - sum_h2 = sum_h2 + hv_obs(k) - enddo - do k=1,nz_data - hv_obs(k)=( sum_h1 / sum_h2 ) * hv_obs(k) - enddo - ! remap increment profile on hv - call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & - nz, hv(1:nz), tmp_val1, h_neglect, h_neglect_edge) - ! add increment to v-velocity on hv - do k=1,nz - v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k) - ! store increment for diagnostics - tmp_v(i,j,k) = tmp_val1(k) - enddo - endif - enddo; enddo + call pass_var(h,G%Domain) ! to ensure reproducibility + + ! add increments to u + hu(:) = 0.0 + tmp_u(:,:,:) = 0.0 ! diagnostics + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) == 1) then + do k=1,nz_data + ! get u increment + tmp_val2(k) = CS%Inc_u%p(i,j,k) + ! get the h and h_obs at u points + hu_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i+1,j,k) ) + enddo + do k=1,nz + hu(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k) ) + enddo + ! account for different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1 + hu(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2 + hu_obs(k) + enddo + do k=1,nz_data + hu_obs(k) = ( sum_h1 / sum_h2 ) * hu_obs(k) + enddo + ! remap increment profile on hu + call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & + nz, hu(1:nz), tmp_val1, h_neglect, h_neglect_edge) + ! add increment to u-velocity on hu + do k=1,nz + u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k) + ! store increment for diagnostics + tmp_u(i,j,k) = tmp_val1(k) + enddo + endif + enddo ; enddo + + ! add increments to v + hv(:) = 0.0 + tmp_v(:,:,:) = 0.0 ! diagnostics + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) == 1) then + ! get v increment + do k=1,nz_data + tmp_val2(k) = CS%Inc_v%p(i,j,k) + ! get the h and h_obs at v points + hv_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i,j+1,k) ) + enddo + do k=1,nz + hv(k) = 0.5 * (h(i,j,k) + h(i,j+1,k) ) + enddo + ! account for different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1 + hv(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2 + hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=( sum_h1 / sum_h2 ) * hv_obs(k) + enddo + ! remap increment profile on hv + call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & + nz, hv(1:nz), tmp_val1, h_neglect, h_neglect_edge) + ! add increment to v-velocity on hv + do k=1,nz + v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k) + ! store increment for diagnostics + tmp_v(i,j,k) = tmp_val1(k) + enddo + endif + enddo ; enddo endif ! uv_inc diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 446cdb8b45..10b46f7e7d 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -275,15 +275,15 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! u points CS%num_col_u = 0 ; if (present(Iresttime_u_in)) then - Iresttime_u(:,:) = Iresttime_u_in(:,:) + Iresttime_u(:,:) = Iresttime_u_in(:,:) else do j=G%jsc,G%jec ; do I=G%iscB,G%iecB Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) enddo ; enddo endif do j=G%jsc,G%jec ; do I=G%iscB,G%iecB - if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) & - CS%num_col_u = CS%num_col_u + 1 + if ((Iresttime_u(I,j) > 0.0) .and. (G%mask2dCu(I,j) > 0.0)) & + CS%num_col_u = CS%num_col_u + 1 enddo ; enddo if (CS%num_col_u > 0) then diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 7371ba7009..58d6e3417a 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -293,7 +293,7 @@ end subroutine calculate_CVMix_conv logical function CVMix_conv_is_used(param_file) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_is_used, & - default=.false., do_not_log = .true.) + default=.false., do_not_log=.true.) end function CVMix_conv_is_used diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index f1ac4c926a..aca84f59e2 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -274,7 +274,7 @@ end subroutine compute_ddiff_coeffs logical function CVMix_ddiff_is_used(param_file) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_is_used, & - default=.false., do_not_log = .true.) + default=.false., do_not_log=.true.) end function CVMix_ddiff_is_used diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index eed99ceb3f..db98e063d8 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -326,9 +326,9 @@ logical function CVMix_shear_is_used(param_file) ! Local variables logical :: LMD94, PP81 call get_param(param_file, mdl, "USE_LMD94", LMD94, & - default=.false., do_not_log = .true.) + default=.false., do_not_log=.true.) call get_param(param_file, mdl, "Use_PP81", PP81, & - default=.false., do_not_log = .true.) + default=.false., do_not_log=.true.) CVMix_shear_is_used = (LMD94 .or. PP81) end function CVMix_shear_is_used diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index c3ee727573..fb218a7d67 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -523,7 +523,7 @@ end subroutine calculate_bkgnd_mixing logical function CVMix_bkgnd_is_used(param_file) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters call get_param(param_file, mdl, "USE_CVMix_BACKGROUND", CVMix_bkgnd_is_used, & - default=.false., do_not_log = .true.) + default=.false., do_not_log=.true.) end function CVMix_bkgnd_is_used diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index cc240dfa95..8f82dce3c9 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3370,7 +3370,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',& 'boundary_forcing_heat_tendency', diag%axesTL, Time, & 'Boundary forcing heat tendency', & - 'W m-2', conversion=US%QRZ_T_to_W_m2, v_extensive = .true.) + 'W m-2', conversion=US%QRZ_T_to_W_m2, v_extensive=.true.) if (CS%id_boundary_forcing_heat_tend > 0) then CS%boundary_forcing_tendency_diag = .true. endif diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index 28a9501d51..2928a0c718 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -288,14 +288,14 @@ subroutine flux_init_OCMIP2_CFC(CS, verbosity) ! These calls obtain the indices for the CFC11 and CFC12 flux coupling. They ! can safely be called multiple times. ind_flux(1) = atmos_ocn_coupler_flux('cfc_11_flux', & - flux_type = 'air_sea_gas_flux', implementation='ocmip2', & + flux_type='air_sea_gas_flux', implementation='ocmip2', & param=(/ 9.36e-07, 9.7561e-06 /), & - ice_restart_file = default_ice_restart_file, & - ocean_restart_file = default_ocean_restart_file, & - caller = "register_OCMIP2_CFC", verbosity=verbosity) + ice_restart_file=default_ice_restart_file, & + ocean_restart_file=default_ocean_restart_file, & + caller="register_OCMIP2_CFC", verbosity=verbosity) ind_flux(2) = atmos_ocn_coupler_flux('cfc_12_flux', & flux_type='air_sea_gas_flux', implementation='ocmip2', & - param = (/ 9.36e-07, 9.7561e-06 /), & + param=(/ 9.36e-07, 9.7561e-06 /), & ice_restart_file=default_ice_restart_file, & ocean_restart_file=default_ocean_restart_file, & caller="register_OCMIP2_CFC", verbosity=verbosity) diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 5e8632f1ca..65947fc64c 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -129,8 +129,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) !Register all the generic tracers used and create the list of them. !This can be called by ALL PE's. No array fields allocated. if (.not. g_registered) then - call generic_tracer_register() - g_registered = .true. + call generic_tracer_register() + g_registered = .true. endif @@ -188,29 +188,29 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) g_tracer=>CS%g_tracer_list do - call g_tracer_get_alias(g_tracer,g_tracer_name) - - call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) - call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname) - call g_tracer_get_values(g_tracer,g_tracer_name,'units',units ) - - !!nnz: MOM field is 3D. Does this affect performance? Need it be override field? - tr_ptr => tr_field(:,:,:,1) - ! Register prognastic tracer for horizontal advection, diffusion, and restarts. - if (g_tracer_is_prog(g_tracer)) then - call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & - name=g_tracer_name, longname=longname, units=units, & - registry_diags=.false., & !### CHANGE TO TRUE? - restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) - else - call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & - restart_CS, longname=longname, units=units) - endif - - !traverse the linked list till hit NULL - call g_tracer_get_next(g_tracer, g_tracer_next) - if (.NOT. associated(g_tracer_next)) exit - g_tracer=>g_tracer_next + call g_tracer_get_alias(g_tracer,g_tracer_name) + + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) + call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname) + call g_tracer_get_values(g_tracer,g_tracer_name,'units',units ) + + !!nnz: MOM field is 3D. Does this affect performance? Need it be override field? + tr_ptr => tr_field(:,:,:,1) + ! Register prognastic tracer for horizontal advection, diffusion, and restarts. + if (g_tracer_is_prog(g_tracer)) then + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + name=g_tracer_name, longname=longname, units=units, & + registry_diags=.false., & !### CHANGE TO TRUE? + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + else + call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & + restart_CS, longname=longname, units=units) + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next enddo @@ -281,70 +281,70 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, if (.not.restart .or. (CS%tracers_may_reinit .and. & .not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then - if (g_tracer%requires_src_info ) then - call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& - "initializing generic tracer "//trim(g_tracer_name)//& - " using MOM_initialize_tracer_from_Z ") - - call MOM_initialize_tracer_from_Z(h, tr_ptr, G, GV, US, param_file, & - src_file = g_tracer%src_file, & - src_var_nam = g_tracer%src_var_name, & - src_var_unit_conversion = g_tracer%src_var_unit_conversion,& - src_var_record = g_tracer%src_var_record, & - src_var_gridspec = g_tracer%src_var_gridspec ) - - !Check/apply the bounds for each g_tracer - do k=1,nk ; do j=jsc,jec ; do i=isc,iec - if (tr_ptr(i,j,k) /= CS%tracer_land_val) then - if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min - !Jasmin does not want to apply the maximum for now - !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max - endif - enddo ; enddo ; enddo - - !jgj: Reset CASED to 0 below K=1 - if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then - do k=2,nk ; do j=jsc,jec ; do i=isc,iec - if (tr_ptr(i,j,k) /= CS%tracer_land_val) then - tr_ptr(i,j,k) = 0.0 - endif - enddo ; enddo ; enddo - endif - elseif(.not. g_tracer%requires_restart) then + if (g_tracer%requires_src_info ) then + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initializing generic tracer "//trim(g_tracer_name)//& + " using MOM_initialize_tracer_from_Z ") + + call MOM_initialize_tracer_from_Z(h, tr_ptr, G, GV, US, param_file, & + src_file = g_tracer%src_file, & + src_var_nam = g_tracer%src_var_name, & + src_var_unit_conversion = g_tracer%src_var_unit_conversion,& + src_var_record = g_tracer%src_var_record, & + src_var_gridspec = g_tracer%src_var_gridspec ) + + !Check/apply the bounds for each g_tracer + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min + !Jasmin does not want to apply the maximum for now + !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max + endif + enddo ; enddo ; enddo + + !jgj: Reset CASED to 0 below K=1 + if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then + do k=2,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + tr_ptr(i,j,k) = 0.0 + endif + enddo ; enddo ; enddo + endif + elseif(.not. g_tracer%requires_restart) then !Do nothing for this tracer, it is initialized by the tracer package call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& "skip initialization of generic tracer "//trim(g_tracer_name)) - else !Do it old way if the tracer is not registered to start from a specific source file. + else !Do it old way if the tracer is not registered to start from a specific source file. !This path should be deprecated if all generic tracers are required to start from specified sources. - if (len_trim(CS%IC_file) > 0) then - ! Read the tracer concentrations from a netcdf file. - if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, & - "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file) - if (CS%Z_IC_file) then - OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US) - if (.not.OK) then - OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US) - if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//& - "Unable to read "//trim(g_tracer_name)//" from "//& - trim(CS%IC_file)//".") + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, & + "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//& + "Unable to read "//trim(g_tracer_name)//" from "//& + trim(CS%IC_file)//".") + endif + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initialized generic tracer "//trim(g_tracer_name)//& + " using Generic Tracer File on Z: "//CS%IC_file) + else + ! native grid + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) + call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain) endif - call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& - "initialized generic tracer "//trim(g_tracer_name)//& - " using Generic Tracer File on Z: "//CS%IC_file) else - ! native grid - call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& - "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//& - " for tracer "//trim(g_tracer_name)) - call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain) + call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//& + "check Generic Tracer IC filename "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) endif - else - call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//& - "check Generic Tracer IC filename "//trim(CS%IC_file)//& - " for tracer "//trim(g_tracer_name)) - endif - endif + endif endif !traverse the linked list till hit NULL @@ -459,21 +459,21 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! g_tracer=>CS%g_tracer_list do - if (_ALLOCATED(g_tracer%trunoff)) then - call g_tracer_get_alias(g_tracer,g_tracer_name) - call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array) - call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array) - call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array) - !nnz: Why is fluxes%river = 0? - runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * & - US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:) - stf_array = stf_array + runoff_tracer_flux_array - endif - - !traverse the linked list till hit NULL - call g_tracer_get_next(g_tracer, g_tracer_next) - if (.NOT. associated(g_tracer_next)) exit - g_tracer=>g_tracer_next + if (_ALLOCATED(g_tracer%trunoff)) then + call g_tracer_get_alias(g_tracer,g_tracer_name) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array) + !nnz: Why is fluxes%river = 0? + runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * & + US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:) + stf_array = stf_array + runoff_tracer_flux_array + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer => g_tracer_next enddo diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 7296640560..8efb341ec3 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -126,7 +126,7 @@ logical function lateral_boundary_diffusion_init(Time, G, GV, param_file, diag, "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ,& - check_reconstruction = .false., check_remapping = .false., answers_2018 = .false.) + check_reconstruction=.false., check_remapping=.false., answers_2018=.false.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) call get_param(param_file, mdl, "LBD_DEBUG", CS%debug, & "If true, write out verbose debugging data in the LBD module.", & @@ -382,9 +382,9 @@ subroutine unique(val, n, val_unique, val_max) max_val = MAXVAL(val) i = 0 do while (min_valmin_val) - tmp(i) = min_val + i = i+1 + min_val = MINVAL(val, mask=val>min_val) + tmp(i) = min_val enddo ii = i if (limit) then @@ -745,8 +745,8 @@ logical function near_boundary_unit_tests( verbose ) allocate(CS) ! fill required fields in CS CS%linear=.false. - call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation = .true. ,& - check_reconstruction = .true., check_remapping = .true.) + call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation=.true. ,& + check_reconstruction=.true., check_remapping=.true.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) CS%H_subroundoff = 1.0E-20 CS%debug=.false. diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 08361ba9af..7f1f39bbc0 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -167,11 +167,11 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", CS%interior_only, & "If true, only applies neutral diffusion in the ocean interior."//& "That is, the algorithm will exclude the surface and bottom"//& - "boundary layers.", default = .false.) + "boundary layers.", default=.false.) call get_param(param_file, mdl, "NDIFF_USE_UNMASKED_TRANSPORT_BUG", CS%use_unmasked_transport_bug, & "If true, use an older form for the accumulation of neutral-diffusion "//& "transports that were unmasked, as used prior to Jan 2018. This is not "//& - "recommended.", default = .false.) + "recommended.", default=.false.) ! Initialize and configure remapping if ( .not.CS%continuous_reconstruction ) then diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index f3d77d2ef3..5012cc4d66 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -1498,7 +1498,7 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic call get_param(param_File, mdl, "RECALC_NEUTRAL_SURF", CS%recalc_neutral_surf, & "If true, then recalculate the neutral surfaces if the \n"//& "diffusive CFL is exceeded. If false, assume that the \n"//& - "positions of the surfaces do not change \n", default = .false.) + "positions of the surfaces do not change \n", default=.false.) CS%ML_KhTR_scale = 1.0 if (CS%Diffuse_ML_interior) then call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, & diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 6214f2d095..64fb31f68d 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -223,7 +223,7 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& "toward some specified surface state with a rate "//& - "given by FLUXCONST.", default= .false.) + "given by FLUXCONST.", default=.false.) if (CS%restorebuoy) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//& diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 8c980edc8e..54985e35fc 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -228,7 +228,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& "toward some specified surface state with a rate "//& - "given by FLUXCONST.", default= .false.) + "given by FLUXCONST.", default=.false.) if (CS%restorebuoy) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//&