diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 93eabff63d..6d480b27ed 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3553,14 +3553,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif endif elseif (CS%use_RK2) then - call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, Time, G, GV, US, & - param_file, diag, CS%dyn_unsplit_RK2_CSp, & + call initialize_dyn_unsplit_RK2(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, & + US, param_file, diag, CS%dyn_unsplit_RK2_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%OBC, & CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, & CS%ntrunc, cont_stencil=CS%cont_stencil) else - call initialize_dyn_unsplit(CS%u, CS%v, CS%h, Time, G, GV, US, & - param_file, diag, CS%dyn_unsplit_CSp, & + call initialize_dyn_unsplit(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, & + US, param_file, diag, CS%dyn_unsplit_CSp, & CS%ADp, CS%CDp, MOM_internal_state, CS%OBC, & CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, CS%visc, dirs, & CS%ntrunc, cont_stencil=CS%cont_stencil) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index c012d15c2f..d11af637a1 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1552,7 +1552,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) - if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS) if (CS%use_tides) then call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) HA_CSp => CS%HA_CSp diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 8a373e61c8..9bfbff5191 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -1435,7 +1435,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) - if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS) if (CS%use_tides) then call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) HA_CSp => CS%HA_CSp diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index bce0c4026a..d4d3356c3d 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -576,7 +576,7 @@ subroutine register_restarts_dyn_unsplit(HI, GV, param_file, CS) end subroutine register_restarts_dyn_unsplit !> Initialize parameters and allocate memory associated with the unsplit dynamics module. -subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS, & +subroutine initialize_dyn_unsplit(u, v, h, tv, Time, G, GV, US, param_file, diag, CS, & Accel_diag, Cont_diag, MIS, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, cont_stencil) @@ -589,6 +589,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type type(time_type), target, intent(in) :: Time !< The current model time. type(param_file_type), intent(in) :: param_file !< A structure to parse !! for run-time parameters. @@ -710,7 +711,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) - if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp) if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, CS%ADp, & CS%SAL_CSp, CS%tides_CSp) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index dd3df7bb3a..7bdae9ab20 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -526,7 +526,7 @@ subroutine register_restarts_dyn_unsplit_RK2(HI, GV, param_file, CS) end subroutine register_restarts_dyn_unsplit_RK2 !> Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module. -subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag, CS, & +subroutine initialize_dyn_unsplit_RK2(u, v, h, tv, Time, G, GV, US, param_file, diag, CS, & Accel_diag, Cont_diag, MIS, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & visc, dirs, ntrunc, cont_stencil) @@ -536,6 +536,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< The zonal velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type type(time_type), target, intent(in) :: Time !< The current model time. type(param_file_type), intent(in) :: param_file !< A structure to parse !! for run-time parameters. @@ -673,7 +674,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv) - if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp) if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, CS%ADp, & CS%SAL_CSp, CS%tides_CSp) diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 5aa822a000..c9e4bc015e 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -3,7 +3,7 @@ module MOM_interface_heights ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol +use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol, int_density_dz use MOM_debugging, only : hchksum use MOM_error_handler, only : MOM_error, FATAL use MOM_EOS, only : calculate_density, average_specific_vol, EOS_type, EOS_domain @@ -20,7 +20,7 @@ module MOM_interface_heights public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple public calc_derived_thermo public convert_MLD_to_ML_thickness -public find_rho_bottom, find_col_avg_SpV +public find_rho_bottom, find_col_avg_SpV, find_col_mass !> Calculates the heights of the free surface or all interfaces from layer thicknesses. interface find_eta @@ -73,7 +73,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. ! dZ_ref is 0 unless the optional argument dZref is present. - integer i, j, k, isv, iev, jsv, jev, nz, halo + integer :: i, j, k, isv, iev, jsv, jev, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -191,7 +191,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. ! dZ_ref is 0 unless the optional argument dZref is present. - integer i, j, k, is, ie, js, je, nz, halo + integer :: i, j, k, is, ie, js, je, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo @@ -345,7 +345,7 @@ subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size) real :: I_rho ! The inverse of the Boussiensq reference density [R-1 ~> m3 kg-1] real :: SpV_lay(SZK_(GV)) ! The inverse of the layer target potential densities [R-1 ~> m3 kg-1] character(len=128) :: mesg ! A string for error messages - integer i, j, k, is, ie, js, je, nz, halo + integer :: i, j, k, is, ie, js, je, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -391,6 +391,71 @@ subroutine find_col_avg_SpV(h, SpV_avg, tv, G, GV, US, halo_size) end subroutine find_col_avg_SpV +!> Calculate the integrated mass of the water column. +subroutine find_col_mass(h, tv, G, GV, US, mass, p_bot, p_surf) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass !< Integrated mass of the water column + !! [R Z ~> kg m-2] + real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: p_bot !< Bottom pressure = g * mass + psurf + !! [R L2 T-2 ~> Pa] + real, dimension(:,:), optional, pointer :: p_surf !< A pointer to surface pressure + !! [R L2 T-2 ~> Pa] + + ! Local variables + real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1] + real, dimension(SZI_(G),SZJ_(G)) :: & + z_top, & ! Height of the top of a layer [Z ~> m]. + z_bot, & ! Height of the bottom of a layer [Z ~> m]. + dp ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa]. + integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isq = G%iscB ; ieq = G%iecB ; jsq = G%jscB ; jeq = G%jecB + nz = GV%ke + + do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo + if (GV%Boussinesq) then + if (associated(tv%eqn_of_state)) then + I_gEarth = 1.0 / GV%g_Earth + do j=jsq,jeq+1 ; do i=isq,ieq+1 ; z_bot(i,j) = 0.0 ; enddo ; enddo + do k=1,nz + ! NOTE: int_density_z expects z_top and z_bot values from [ij]sq to [ij]eq+1 + do j=jsq,jeq+1 ; do i=isq,ieq+1 + z_top(i,j) = z_bot(i,j) + z_bot(i,j) = z_top(i,j) - GV%H_to_Z * h(i,j,k) + enddo ; enddo + call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, & + G%HI, tv%eqn_of_state, US, dp) + do j=js,je ; do i=is,ie + mass(i,j) = mass(i,j) + dp(i,j) * I_gEarth + enddo ; enddo + enddo + else + do k=1,nz ; do j=js,je ; do i=is,ie + mass(i,j) = mass(i,j) + (GV%H_to_Z * GV%Rlay(k)) * h(i,j,k) + enddo ; enddo ; enddo + endif + else + do k=1,nz ; do j=js,je ; do i=is,ie + mass(i,j) = mass(i,j) + GV%H_to_RZ * h(i,j,k) + enddo ; enddo ; enddo + endif + + if (present(p_bot)) then + do j=js,je ; do i=is,ie + p_bot(i,j) = GV%g_Earth * mass(i,j) + enddo ; enddo + if (present(p_surf) .and. associated(p_surf)) then ; do j=js,je ; do i=is,ie + p_bot(i,j) = p_bot(i,j) + p_surf(i,j) + enddo ; enddo ; endif + endif + +end subroutine find_col_mass !> Determine the in situ density averaged over a specified distance from the bottom, !! calculating it as the inverse of the mass-weighted average specific volume. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index d164363ec4..6c220c79cf 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -24,7 +24,7 @@ module MOM_diagnostics use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_interface_heights, only : find_eta +use MOM_interface_heights, only : find_eta, find_col_mass use MOM_spatial_means, only : global_area_mean, global_layer_mean use MOM_spatial_means, only : global_volume_mean, global_area_integral use MOM_tracer_registry, only : tracer_registry_type, post_tracer_transport_diagnostics @@ -39,7 +39,6 @@ module MOM_diagnostics #include public calculate_diagnostic_fields, register_time_deriv, write_static_fields -public find_eta public register_surface_diags, post_surface_dyn_diags, post_surface_thermo_diags public register_transport_diags, post_transport_diagnostics public MOM_diagnostics_init, MOM_diagnostics_end @@ -903,7 +902,6 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'. ! This is the column mass multiplied by gravity plus the pressure ! at the ocean surface [R L2 T-2 ~> Pa]. - dpress, & ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa]. tr_int ! vertical integral of a tracer times density, ! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2]. real :: IG_Earth ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1]. @@ -943,52 +941,14 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) call post_data(CS%id_col_ht, z_bot, CS%diag) endif - ! NOTE: int_density_z expects z_top and z_btm values from [ij]sq to [ij]eq+1 if (CS%id_col_mass > 0 .or. CS%id_pbo > 0) then - do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo - if (GV%Boussinesq) then - if (associated(tv%eqn_of_state)) then - IG_Earth = 1.0 / GV%g_Earth - do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1 - z_bot(i,j) = 0.0 - enddo ; enddo - do k=1,nz - do j=G%jscB,G%jecB+1 ; do i=G%iscB,G%iecB+1 - z_top(i,j) = z_bot(i,j) - z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k) - enddo ; enddo - call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, & - G%HI, tv%eqn_of_state, US, dpress) - do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth - enddo ; enddo - enddo - else - do k=1,nz ; do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + (GV%H_to_Z*GV%Rlay(k))*h(i,j,k) - enddo ; enddo ; enddo - endif - else - do k=1,nz ; do j=js,je ; do i=is,ie - mass(i,j) = mass(i,j) + GV%H_to_RZ*h(i,j,k) - enddo ; enddo ; enddo - endif - if (CS%id_col_mass > 0) then - call post_data(CS%id_col_mass, mass, CS%diag) - endif if (CS%id_pbo > 0) then - do j=js,je ; do i=is,ie ; btm_pres(i,j) = 0.0 ; enddo ; enddo - ! 'pbo' is defined as the sea water pressure at the sea floor - ! pbo = (mass * g) + p_surf - ! where p_surf is the sea water pressure at sea water surface. - do j=js,je ; do i=is,ie - btm_pres(i,j) = GV%g_Earth * mass(i,j) - if (associated(p_surf)) then - btm_pres(i,j) = btm_pres(i,j) + p_surf(i,j) - endif - enddo ; enddo + call find_col_mass(h, tv, G, GV, US, mass, btm_pres, p_surf) call post_data(CS%id_pbo, btm_pres, CS%diag) + else + call find_col_mass(h, tv, G, GV, US, mass) endif + if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag) endif end subroutine calculate_vertical_integrals diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index 5b2ba9bad1..045027f05c 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -1,18 +1,22 @@ module MOM_self_attr_load -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE -use MOM_domains, only : pass_var -use MOM_error_handler, only : MOM_error, FATAL, WARNING -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_grid, only : ocean_grid_type -use MOM_io, only : slasher, MOM_read_data -use MOM_load_love_numbers, only : Love_Data -use MOM_obsolete_params, only : obsolete_logical, obsolete_int +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE +use MOM_domains, only : pass_var +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : find_col_mass +use MOM_io, only : MOM_infra_file, MOM_field, vardesc, slasher +use MOM_io, only : create_MOM_file, MOM_read_data, MOM_write_field, var_desc +use MOM_load_love_numbers, only : Love_Data +use MOM_restart, only : is_new_run, MOM_restart_CS use MOM_spherical_harmonics, only : spherical_harmonics_init, spherical_harmonics_end use MOM_spherical_harmonics, only : spherical_harmonics_forward, spherical_harmonics_inverse use MOM_spherical_harmonics, only : sht_CS, order2index, calc_lmax -use MOM_unit_scaling, only : unit_scale_type -use MOM_verticalGrid, only : verticalGrid_type +use MOM_string_functions, only : lowercase +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -39,8 +43,8 @@ module MOM_self_attr_load !< Spherical harmonic transforms (SHT) control structure integer :: sal_sht_Nd !< Maximum degree for spherical harmonic transforms [nondim] - real, allocatable :: ebot_ref(:,:) - !< Reference bottom pressure scaled by Rho_0 and G_Earth[Z ~> m] + real, allocatable :: pbot_ref(:,:) + !< Reference bottom pressure [R L2 T-2 ~> Pa] real, allocatable :: Love_scaling(:) !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers !! [nondim] or [Z T2 L-2 R-1 ~> m Pa-1], depending on the value of use_ppa. @@ -79,7 +83,7 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB if (CS%use_bpa) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - bpa(i,j) = eta(i,j) - CS%ebot_ref(i,j) + bpa(i,j) = eta(i,j) - CS%pbot_ref(i,j) enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 bpa(i,j) = eta(i,j) enddo ; enddo ; endif @@ -176,20 +180,30 @@ subroutine calc_love_scaling(rhoW, rhoE, grav, CS) end subroutine calc_love_scaling !> This subroutine initializes the self-attraction and loading control structure. -subroutine SAL_init(G, GV, US, param_file, CS) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. - type(SAL_CS), intent(inout) :: CS !< Self-attraction and loading control structure - +subroutine SAL_init(h, tv, G, GV, US, param_file, CS, restart_CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(SAL_CS), intent(inout) :: CS !< Self-attraction and loading control structure + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(MOM_restart_CS), optional, intent(in) :: restart_CS !< MOM restart control structure ! Local variables # include "version_variable.h" character(len=40) :: mdl = "MOM_self_attr_load" ! This module's name. integer :: lmax ! Total modes of the real spherical harmonics [nondim] real :: rhoE ! The average density of Earth [R ~> kg m-3]. - character(len=200) :: filename, ebot_ref_file, inputdir ! Strings for file/path - character(len=200) :: ebot_ref_varname ! Variable name in file + character(len=20) :: bpa_config ! String for reference bottom pressure config option + real :: tmp(G%isd:G%ied, G%jsd:G%jed) ! Temporary field storing mass returned by find_col_mass + ! [R Z ~> kg m-2] + logical :: restart_sim ! If true, this is a restart run + character(len=200) :: filename, ref_pbot_file, inputdir ! Strings for file/path + character(len=200) :: ref_pbot_varname ! Variable name in file + type(MOM_infra_file) :: IO_handle ! used to write ref_pbot file + type(vardesc) :: vars(1) ! used to write ref_pbot file + type(MOM_field) :: fields(1) ! used to write ref_pbot file logical :: calculate_sal, tides, use_tidal_sal_file integer :: tides_answer_date ! Recover old answers with tides real :: sal_scalar_value ! Scaling SAL factors [nondim] @@ -204,40 +218,62 @@ subroutine SAL_init(G, GV, US, param_file, CS) call get_param(param_file, '', "CALCULATE_SAL", calculate_sal, default=tides, do_not_log=.True.) if (.not. calculate_sal) return - if (tides) then - call get_param(param_file, '', "USE_PREVIOUS_TIDES", CS%use_tidal_sal_prev, & - default=.false., do_not_log=.True.) - call get_param(param_file, '', "TIDAL_SAL_FROM_FILE", use_tidal_sal_file, & - default=.false., do_not_log=.True.) - call get_param(param_file, '', "TIDES_ANSWER_DATE", tides_answer_date, & - default=20230630, do_not_log=.True.) - endif - call get_param(param_file, mdl, "SAL_USE_BPA", CS%use_bpa, & "If true, use bottom pressure anomaly to calculate self-attraction and "// & "loading (SAL). Otherwise sea surface height anomaly is used, which is "// & - "only correct for homogenous flow.", default=.False.) + "only accurate for uniform density fluid.", default=.False.) if (CS%use_bpa) then + allocate(CS%pbot_ref(isd:ied, jsd:jed), source=0.0) + call get_param(param_file, mdl, "SAL_REF_PBOT_CONFIG", bpa_config, default="file", & + do_not_log=.True.) + restart_sim = .False. ; if (present(restart_CS)) restart_sim = (.not. is_new_run(restart_CS)) + if (restart_sim .and. (trim(lowercase(bpa_config))/='file')) then + call MOM_error(WARNING, "SAL_init: 'file' is not used by SAL_PBOT_REF_CONFIG for a restart "//& + "run, SAL_PBOT_REF_CONFIG is reset to 'file'.") + bpa_config = 'file' + endif + call get_param(param_file, mdl, "SAL_REF_PBOT_CONFIG", bpa_config, & + "A string that determines how the reference bottom pressure for SAL "//& + "is specified:\n"//& + "\t init - calculated by thickness, temperature and salinity from \n"//& + "\t initialization and assuming surface pressure is zero.\n"//& + "\t This option can only be used by new simulations.\n"//& + "\t file - read from the file specified by REF_PBOT_FILE.", & + default="file", do_not_read=.True.) call get_param(param_file, '', "INPUTDIR", inputdir, default=".", do_not_log=.True.) - inputdir = slasher(inputdir) - call get_param(param_file, mdl, "REF_BOT_PRES_FILE", ebot_ref_file, & + call get_param(param_file, mdl, "REF_PBOT_FILE", ref_pbot_file, & "Reference bottom pressure file used by self-attraction and loading (SAL).", & default="pbot.nc") - call get_param(param_file, mdl, "REF_BOT_PRES_VARNAME", ebot_ref_varname, & - "The name of the variable in REF_BOT_PRES_FILE with reference bottom "//& - "pressure. The variable should have the unit of Pa.", & - default="pbot") - filename = trim(inputdir)//trim(ebot_ref_file) - call log_param(param_file, mdl, "INPUTDIR/REF_BOT_PRES_FILE", filename) - - allocate(CS%ebot_ref(isd:ied, jsd:jed), source=0.0) - call MOM_read_data(filename, trim(ebot_ref_varname), CS%ebot_ref, G%Domain,& - scale=US%Pa_to_RL2_T2) - call pass_var(CS%ebot_ref, G%Domain) + call get_param(param_file, mdl, "REF_PBOT_VARNAME", ref_pbot_varname, & + "The name of the variable in REF_PBOT_FILE with reference bottom "//& + "pressure. The variable should have the unit of Pa.", default="pbot") + filename = trim(slasher(inputdir))//trim(ref_pbot_file) + call log_param(param_file, mdl, "INPUTDIR/REF_PBOT_FILE", filename) + select case (trim(lowercase(bpa_config))) + case ("file") + call MOM_read_data(filename, trim(ref_pbot_varname), CS%pbot_ref, G%Domain,& + scale=US%Pa_to_RL2_T2) + case ("init") + call find_col_mass(h, tv, G, GV, US, tmp, CS%pbot_ref) + ! Write reference bottom pressure file + vars(1) = var_desc(trim(ref_pbot_varname), units="Pa", & + longname="Reference bottom pressure", & + hor_grid='h', z_grid='1', t_grid='1') + call create_MOM_file(IO_handle, trim(filename), vars, 1, fields, G=G) + call MOM_write_field(IO_handle, fields(1), G%Domain, CS%pbot_ref, unscale=US%RL2_T2_to_Pa) + call IO_handle%close() + case default + call MOM_error(FATAL, "SAL_init: Unsupported SAL_PBOT_REF_CONFIG option "//trim(bpa_config)) + end select + call pass_var(CS%pbot_ref, G%Domain) endif + call get_param(param_file, '', "TIDES_ANSWER_DATE", tides_answer_date, default=20230630, & + do_not_log=.True.) ! used to check SAL_USE_BPA if (tides_answer_date<=20250131 .and. CS%use_bpa) & call MOM_error(FATAL, trim(mdl) // ", SAL_init: SAL_USE_BPA needs to be false to recover "//& "tide answers before 20250131.") + call get_param(param_file, '', "TIDAL_SAL_FROM_FILE", use_tidal_sal_file, default=.false., & + do_not_log=.True.) ! used to set default of SAL_SCALAR_APPROX call get_param(param_file, mdl, "SAL_SCALAR_APPROX", CS%use_sal_scalar, & "If true, use the scalar approximation to calculate self-attraction and "//& "loading.", default=tides .and. (.not. use_tidal_sal_file)) @@ -250,6 +286,8 @@ subroutine SAL_init(G, GV, US, param_file, CS) "SAL_SCALAR_APPROX is true or USE_PREVIOUS_TIDES is true.", default=0.0, & units="m m-1", do_not_log=.not.(CS%use_sal_scalar .or. CS%use_tidal_sal_prev), & old_name='TIDE_SAL_SCALAR_VALUE') + call get_param(param_file, '', "USE_PREVIOUS_TIDES", CS%use_tidal_sal_prev, & + default=.false., do_not_log=.True.) call get_param(param_file, mdl, "SAL_HARMONICS", CS%use_sal_sht, & "If true, use the online spherical harmonics method to calculate "//& "self-attraction and loading.", default=.false.) @@ -300,7 +338,7 @@ subroutine SAL_end(CS) type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call !! to SAL_init; it is deallocated here. - if (allocated(CS%ebot_ref)) deallocate(CS%ebot_ref) + if (allocated(CS%pbot_ref)) deallocate(CS%pbot_ref) if (CS%use_sal_sht) then if (allocated(CS%Love_scaling)) deallocate(CS%Love_scaling)