diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 03d23ebf7b..ebe074f0b4 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -877,25 +877,28 @@ end subroutine ocean_model_init_sfc !WGA subroutine ocean_model_flux_init(OS) - type(ocean_state_type), pointer :: OS + type(ocean_state_type), optional, pointer :: OS + integer :: dummy character(len=128) :: default_ice_restart_file, default_ocean_restart_file character(len=40) :: mod = "ocean_model_flux_init" ! This module's name. - type(param_file_type) :: param_file !< A structure to parse for run-time parameters type(directories) :: dirs_tmp ! A structure containing several relevant directory paths. logical :: use_OCMIP_CFCs, use_MOM_generic_tracer + logical :: OS_is_set + + OS_is_set = .false. ; if (present(OS)) OS_is_set = associated(OS) call get_MOM_Input(param_file, dirs_tmp, check_params=.false.) call get_param(param_file, mod, "USE_OCMIP2_CFC", use_OCMIP_CFCs, & - default=.false.) + default=.false., do_not_log=.true.) call get_param(param_file, mod, "USE_generic_tracer", use_MOM_generic_tracer,& - default=.false.) + default=.false., do_not_log=.true.) call close_param_file(param_file, quiet_close=.true.) - if(.not.associated(OS)) then + if(.not.OS_is_set) then if (use_OCMIP_CFCs)then default_ice_restart_file = 'ice_ocmip2_cfc.res.nc' default_ocean_restart_file = 'ocmip2_cfc.res.nc' @@ -917,7 +920,7 @@ subroutine ocean_model_flux_init(OS) if (use_MOM_generic_tracer) then #ifdef _USE_GENERIC_TRACER - call MOM_generic_flux_init + call MOM_generic_flux_init() #else call MOM_error(FATAL, & "call_tracer_register: use_MOM_generic_tracer=.true. BUT not compiled with _USE_GENERIC_TRACER") diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index f522d07ae9..a0d56b58a0 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -154,6 +154,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & integer :: regrid_iterations logical :: Analytic_FV_PGF, obsol_test logical :: convert + logical :: just_read ! If true, only read the parameters because this + ! is a run from a restart file; this option + ! allows the use of Fatal unused parameters. type(EOS_type), pointer :: eos => NULL() logical :: debug ! indicates whether to write debugging output ! This include declares and sets the variable "version". @@ -174,6 +177,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. + just_read = .not.new_sim + call get_param(PF, mod, "INPUTDIR", inputdir, & "The directory in which input files are found.", default=".") inputdir = slasher(inputdir) @@ -189,225 +194,264 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & !==================================================================== if (new_sim) then -! This block initializes all of the fields internally. ! call MOM_mesg("Run initialized internally.", 3) if (present(Time_in)) Time = Time_in ! Otherwise leave Time at its input value. - call get_param(PF, mod, "INIT_LAYERS_FROM_Z_FILE", from_Z_file, & - "If true, intialize the layer thicknesses, temperatures, \n"//& - "and salnities from a Z-space file on a latitude- \n"//& - "longitude grid.", default=.false.) ! h will be converted from m to H below h(:,:,:) = GV%Angstrom_z + endif + + ! The remaining initialization calls are done, regardless of whether the + ! fields are actually initialized here (if just_read=.false.) or whether it + ! is just to make sure that all valid parameters are read to enable the + ! detection of unused parameters. + call get_param(PF, mod, "INIT_LAYERS_FROM_Z_FILE", from_Z_file, & + "If true, intialize the layer thicknesses, temperatures, \n"//& + "and salnities from a Z-space file on a latitude- \n"//& + "longitude grid.", default=.false., do_not_log=just_read) - if (from_Z_file) then + if (from_Z_file) then ! Initialize thickness and T/S from z-coordinate data in a file. - if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//& - "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true") + if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//& + "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true") - call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) - call pass_var(h, G%Domain) + call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, just_read_params=just_read) - else + else ! Initialize thickness, h. - call get_param(PF, mod, "THICKNESS_CONFIG", config, & - "A string that determines how the initial layer \n"//& - "thicknesses are specified for a new run: \n"//& - " \t file - read interface heights from the file specified \n"//& - " \t thickness_file - read thicknesses from the file specified \n"//& - " \t\t by (THICKNESS_FILE).\n"//& - " \t coord - determined by ALE coordinate.\n"//& - " \t uniform - uniform thickness layers evenly distributed \n"//& - " \t\t between the surface and MAXIMUM_DEPTH. \n"//& - " \t DOME - use a slope and channel configuration for the \n"//& - " \t\t DOME sill-overflow test case. \n"//& - " \t ISOMIP - use a configuration for the \n"//& - " \t\t ISOMIP test case. \n"//& - " \t benchmark - use the benchmark test case thicknesses. \n"//& - " \t search - search a density profile for the interface \n"//& - " \t\t densities. This is not yet implemented. \n"//& - " \t circle_obcs - the circle_obcs test case is used. \n"//& - " \t DOME2D - 2D version of DOME initialization. \n"//& - " \t adjustment2d - TBD AJA. \n"//& - " \t sloshing - TBD AJA. \n"//& - " \t seamount - TBD AJA. \n"//& - " \t soliton - Equatorial Rossby soliton. \n"//& - " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& - " \t USER - call a user modified routine.", & - fail_if_missing=.true.) - select case (trim(config)) - case ("file"); call initialize_thickness_from_file(h, G, GV, PF, .false.) - case ("thickness_file"); call initialize_thickness_from_file(h, G, GV, PF, .true.) - case ("coord") - if (useALE) then - call ALE_initThicknessToCoord( ALE_CSp, G, GV, h ) - else - call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//& - "for THICKNESS_CONFIG of 'coord'") - endif - case ("uniform"); call initialize_thickness_uniform(h, G, GV, PF) - case ("DOME"); call DOME_initialize_thickness(h, G, GV, PF) - case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, PF, tv) - case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, PF, & - tv%eqn_of_state, tv%P_Ref) - case ("search"); call initialize_thickness_search - case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF) - case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, PF) - case ("external_gwave"); call external_gwave_initialize_thickness(h, G, PF) - case ("DOME2D"); call DOME2d_initialize_thickness(h, G, GV, PF) - case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, PF) - case ("sloshing"); call sloshing_initialize_thickness(h, G, GV, PF) - case ("seamount"); call seamount_initialize_thickness(h, G, GV, PF) - case ("soliton"); call soliton_initialize_thickness(h, G) - case ("phillips"); call Phillips_initialize_thickness(h, G, GV, PF) - case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, PF) - case ("USER"); call user_initialize_thickness(h, G, PF, tv%T) - case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& - "Unrecognized layer thickness configuration "//trim(config)) - end select - call pass_var(h, G%Domain) - -! Initialize temperature and salinity (T and S). - if ( use_temperature ) then - call get_param(PF, mod, "TS_CONFIG", config, & - "A string that determines how the initial tempertures \n"//& - "and salinities are specified for a new run: \n"//& - " \t file - read velocities from the file specified \n"//& - " \t\t by (TS_FILE). \n"//& - " \t fit - find the temperatures that are consistent with \n"//& - " \t\t the layer densities and salinity S_REF. \n"//& - " \t TS_profile - use temperature and salinity profiles \n"//& - " \t\t (read from TS_FILE) to set layer densities. \n"//& - " \t benchmark - use the benchmark test case T & S. \n"//& - " \t linear - linear in logical layer space. \n"//& - " \t DOME2D - 2D DOME initialization. \n"//& - " \t ISOMIP - ISOMIP initialization. \n"//& - " \t adjustment2d - TBD AJA. \n"//& - " \t sloshing - TBD AJA. \n"//& - " \t seamount - TBD AJA. \n"//& - " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& - " \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//& - " \t SCM_CVmix_tests - used in the SCM CVmix tests.\n"//& - " \t USER - call a user modified routine.", & - fail_if_missing=.true.) -! " \t baroclinic_zone - an analytic baroclinic zone. \n"//& - select case (trim(config)) - case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, GV, PF, eos, tv%P_Ref) - case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, G, PF) - case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, & - G, GV, PF, eos, tv%P_Ref) - case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, G, PF) - case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, G, PF) - case ("DOME2D"); call DOME2d_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, PF, eos) - case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, GV, PF, eos) - case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, PF, eos) - case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, & - tv%S, h, G, PF) - case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, & - tv%S, h, G, PF, eos) - case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, & - tv%S, h, G, GV, PF, eos) - case ("rossby_front"); call Rossby_front_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, PF, eos) - case ("SCM_ideal_hurr"); call SCM_idealized_hurricane_TS_init ( tv%T, & - tv%S, h, G, GV, PF) - case ("SCM_CVmix_tests"); call SCM_CVmix_tests_TS_init (tv%T, & - tv%S, h, G, GV, PF) - case ("dense"); call dense_water_initialize_TS(G, GV, PF, eos, tv%T, tv%S, h) - case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, PF, eos) - case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& - "Unrecognized Temp & salt configuration "//trim(config)) - end select - endif - endif ! not from_Z_file. - -! Initialize velocity components, u and v - call get_param(PF, mod, "VELOCITY_CONFIG", config, & - "A string that determines how the initial velocities \n"//& - "are specified for a new run: \n"//& - " \t file - read velocities from the file specified \n"//& - " \t\t by (VELOCITY_FILE). \n"//& - " \t zero - the fluid is initially at rest. \n"//& - " \t uniform - the flow is uniform (determined by\n"//& - " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//& - " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& - " \t soliton - Equatorial Rossby soliton.\n"//& - " \t USER - call a user modified routine.", default="zero") + call get_param(PF, mod, "THICKNESS_CONFIG", config, & + "A string that determines how the initial layer \n"//& + "thicknesses are specified for a new run: \n"//& + " \t file - read interface heights from the file specified \n"//& + " \t thickness_file - read thicknesses from the file specified \n"//& + " \t\t by (THICKNESS_FILE).\n"//& + " \t coord - determined by ALE coordinate.\n"//& + " \t uniform - uniform thickness layers evenly distributed \n"//& + " \t\t between the surface and MAXIMUM_DEPTH. \n"//& + " \t DOME - use a slope and channel configuration for the \n"//& + " \t\t DOME sill-overflow test case. \n"//& + " \t ISOMIP - use a configuration for the \n"//& + " \t\t ISOMIP test case. \n"//& + " \t benchmark - use the benchmark test case thicknesses. \n"//& + " \t search - search a density profile for the interface \n"//& + " \t\t densities. This is not yet implemented. \n"//& + " \t circle_obcs - the circle_obcs test case is used. \n"//& + " \t DOME2D - 2D version of DOME initialization. \n"//& + " \t adjustment2d - TBD AJA. \n"//& + " \t sloshing - TBD AJA. \n"//& + " \t seamount - TBD AJA. \n"//& + " \t soliton - Equatorial Rossby soliton. \n"//& + " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& + " \t USER - call a user modified routine.", & + fail_if_missing=new_sim, do_not_log=just_read) select case (trim(config)) - case ("file"); call initialize_velocity_from_file(u, v, G, PF) - case ("zero"); call initialize_velocity_zero(u, v, G, PF) - case ("uniform"); call initialize_velocity_uniform(u, v, G, PF) - case ("circular"); call initialize_velocity_circular(u, v, G, PF) - case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, PF) - case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, G, GV, PF) - case ("soliton"); call soliton_initialize_velocity(u, v, h, G) - case ("USER"); call user_initialize_velocity(u, v, G, PF) + case ("file"); call initialize_thickness_from_file(h, G, GV, PF, .false., just_read_params=just_read) + case ("thickness_file"); call initialize_thickness_from_file(h, G, GV, PF, .true., just_read_params=just_read) + case ("coord") + if (new_sim .and. useALE) then + call ALE_initThicknessToCoord( ALE_CSp, G, GV, h ) + elseif (new_sim) then + call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//& + "for THICKNESS_CONFIG of 'coord'") + endif + case ("uniform"); call initialize_thickness_uniform(h, G, GV, PF, & + just_read_params=just_read) + case ("DOME"); call DOME_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, PF, tv, & + just_read_params=just_read) + case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, PF, & + tv%eqn_of_state, tv%P_Ref, just_read_params=just_read) + case ("search"); call initialize_thickness_search + case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, & + PF, just_read_params=just_read) + case ("external_gwave"); call external_gwave_initialize_thickness(h, G, & + PF, just_read_params=just_read) + case ("DOME2D"); call DOME2d_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, & + PF, just_read_params=just_read) + case ("sloshing"); call sloshing_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("seamount"); call seamount_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("soliton"); call soliton_initialize_thickness(h, G) + case ("phillips"); call Phillips_initialize_thickness(h, G, GV, PF, & + just_read_params=just_read) + case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, & + PF, just_read_params=just_read) + case ("USER"); call user_initialize_thickness(h, G, PF, tv%T, & + just_read_params=just_read) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& - "Unrecognized velocity configuration "//trim(config)) + "Unrecognized layer thickness configuration "//trim(config)) end select - call pass_vector(u, v, G%Domain) - if (debug) then - call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1) +! Initialize temperature and salinity (T and S). + if ( use_temperature ) then + call get_param(PF, mod, "TS_CONFIG", config, & + "A string that determines how the initial tempertures \n"//& + "and salinities are specified for a new run: \n"//& + " \t file - read velocities from the file specified \n"//& + " \t\t by (TS_FILE). \n"//& + " \t fit - find the temperatures that are consistent with \n"//& + " \t\t the layer densities and salinity S_REF. \n"//& + " \t TS_profile - use temperature and salinity profiles \n"//& + " \t\t (read from TS_FILE) to set layer densities. \n"//& + " \t benchmark - use the benchmark test case T & S. \n"//& + " \t linear - linear in logical layer space. \n"//& + " \t DOME2D - 2D DOME initialization. \n"//& + " \t ISOMIP - ISOMIP initialization. \n"//& + " \t adjustment2d - TBD AJA. \n"//& + " \t sloshing - TBD AJA. \n"//& + " \t seamount - TBD AJA. \n"//& + " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& + " \t SCM_ideal_hurr - used in the SCM idealized hurricane test.\n"//& + " \t SCM_CVmix_tests - used in the SCM CVmix tests.\n"//& + " \t USER - call a user modified routine.", & + fail_if_missing=new_sim, do_not_log=just_read) +! " \t baroclinic_zone - an analytic baroclinic zone. \n"//& + select case (trim(config)) + case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, GV, PF, & + eos, tv%P_Ref, just_read_params=just_read) + case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, G, & + PF, just_read_params=just_read) + case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, & + G, GV, PF, eos, tv%P_Ref, just_read_params=just_read) + case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, & + G, PF, just_read_params=just_read) + case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, G, PF, & + just_read_params=just_read) + case ("DOME2D"); call DOME2d_initialize_temperature_salinity ( tv%T, & + tv%S, h, G, PF, eos, just_read_params=just_read) + case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity ( tv%T, & + tv%S, h, G, GV, PF, eos, just_read_params=just_read) + case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, & + tv%S, h, G, PF, eos, just_read_params=just_read) + case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, & + tv%S, h, G, PF, just_read_params=just_read) + case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, & + tv%S, h, G, PF, eos, just_read_params=just_read) + case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, & + tv%S, h, G, GV, PF, eos, just_read_params=just_read) + case ("rossby_front"); call Rossby_front_initialize_temperature_salinity ( tv%T, & + tv%S, h, G, PF, eos, just_read_params=just_read) + case ("SCM_ideal_hurr"); call SCM_idealized_hurricane_TS_init ( tv%T, & + tv%S, h, G, GV, PF, just_read_params=just_read) + case ("SCM_CVmix_tests"); call SCM_CVmix_tests_TS_init (tv%T, & + tv%S, h, G, GV, PF, just_read_params=just_read) + case ("dense"); call dense_water_initialize_TS(G, GV, PF, eos, tv%T, tv%S, & + h, just_read_params=just_read) + case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, PF, eos, & + just_read_params=just_read) + case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& + "Unrecognized Temp & salt configuration "//trim(config)) + end select endif + endif ! not from_Z_file. + + ! The thicknesses in halo points might be needed to initialize the velocities. + if (new_sim) call pass_var(h, G%Domain) + +! Initialize velocity components, u and v + call get_param(PF, mod, "VELOCITY_CONFIG", config, & + "A string that determines how the initial velocities \n"//& + "are specified for a new run: \n"//& + " \t file - read velocities from the file specified \n"//& + " \t\t by (VELOCITY_FILE). \n"//& + " \t zero - the fluid is initially at rest. \n"//& + " \t uniform - the flow is uniform (determined by\n"//& + " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//& + " \t rossby_front - a mixed layer front in thermal wind balance.\n"//& + " \t soliton - Equatorial Rossby soliton.\n"//& + " \t USER - call a user modified routine.", default="zero", & + do_not_log=just_read) + select case (trim(config)) + case ("file"); call initialize_velocity_from_file(u, v, G, PF, & + just_read_params=just_read) + case ("zero"); call initialize_velocity_zero(u, v, G, PF, & + just_read_params=just_read) + case ("uniform"); call initialize_velocity_uniform(u, v, G, PF, & + just_read_params=just_read) + case ("circular"); call initialize_velocity_circular(u, v, G, PF, & + just_read_params=just_read) + case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, PF, & + just_read_params=just_read) + case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, & + G, GV, PF, just_read_params=just_read) + case ("soliton"); call soliton_initialize_velocity(u, v, h, G) + case ("USER"); call user_initialize_velocity(u, v, G, PF, & + just_read_params=just_read) + case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& + "Unrecognized velocity configuration "//trim(config)) + end select + + if (new_sim) call pass_vector(u, v, G%Domain) + if (debug .and. new_sim) then + call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1) + endif ! Optionally convert the thicknesses from m to kg m-2. This is particularly ! useful in a non-Boussinesq model. - call get_param(PF, mod, "CONVERT_THICKNESS_UNITS", convert, & - "If true, convert the thickness initial conditions from \n"//& - "units of m to kg m-2 or vice versa, depending on whether \n"//& - "BOUSSINESQ is defined. This does not apply if a restart \n"//& - "file is read.", default=.false.) + call get_param(PF, mod, "CONVERT_THICKNESS_UNITS", convert, & + "If true, convert the thickness initial conditions from \n"//& + "units of m to kg m-2 or vice versa, depending on whether \n"//& + "BOUSSINESQ is defined. This does not apply if a restart \n"//& + "file is read.", default=.false., do_not_log=just_read) + if (new_sim) then if (convert .and. .not. GV%Boussinesq) then ! Convert h from m to kg m-2 then to thickness units (H) - call convert_thickness(h, G, GV, PF, tv) + call convert_thickness(h, G, GV, tv) elseif (GV%Boussinesq) then ! Convert h from m to thickness units (H) h(:,:,:) = h(:,:,:)*GV%m_to_H else h(:,:,:) = h(:,:,:)*GV%kg_m2_to_H endif + endif ! Remove the mass that would be displaced by an ice shelf or inverse barometer. - call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, & - "If true, depress the initial surface to avoid huge \n"//& - "tsunamis when a large surface pressure is applied.", & - default=.false.) - call get_param(PF, mod, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & - "If true, cuts way the top of the column for initial conditions\n"//& - "at the depth where the hydrostatic presure matches the imposed\n"//& - "surface pressure which is read from file.", default=.false.) - if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//& - "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True") - if (depress_sfc) call depress_surface(h, G, GV, PF, tv) - if (trim_ic_for_p_surf) call trim_for_ice(PF, G, GV, ALE_CSp, tv, h) - - ! Perhaps we want to run the regridding coordinate generator for multiple - ! iterations here so the initial grid is consistent with the coordinate - if (useALE) then - call get_param(PF, mod, "REGRID_ACCELERATE_INIT", regrid_accelerate, & - "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding\n"//& - "algorithm to push the initial grid to be consistent with the initial\n"//& - "condition. Useful only for state-based and iterative coordinates.", & - default=.false.) - if (regrid_accelerate) then - call get_param(PF, mod, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, & - "The number of regridding iterations to perform to generate\n"//& - "an initial grid that is consistent with the initial conditions.", & - default=1) - + call get_param(PF, mod, "DEPRESS_INITIAL_SURFACE", depress_sfc, & + "If true, depress the initial surface to avoid huge \n"//& + "tsunamis when a large surface pressure is applied.", & + default=.false., do_not_log=just_read) + call get_param(PF, mod, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & + "If true, cuts way the top of the column for initial conditions\n"//& + "at the depth where the hydrostatic presure matches the imposed\n"//& + "surface pressure which is read from file.", default=.false., & + do_not_log=just_read) + if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//& + "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True") + if (depress_sfc) call depress_surface(h, G, GV, PF, tv, just_read_params=just_read) + if (trim_ic_for_p_surf) call trim_for_ice(PF, G, GV, ALE_CSp, tv, h, just_read_params=just_read) + + ! Perhaps we want to run the regridding coordinate generator for multiple + ! iterations here so the initial grid is consistent with the coordinate + if (useALE) then + call get_param(PF, mod, "REGRID_ACCELERATE_INIT", regrid_accelerate, & + "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding\n"//& + "algorithm to push the initial grid to be consistent with the initial\n"//& + "condition. Useful only for state-based and iterative coordinates.", & + default=.false., do_not_log=just_read) + if (regrid_accelerate) then + call get_param(PF, mod, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, & + "The number of regridding iterations to perform to generate\n"//& + "an initial grid that is consistent with the initial conditions.", & + default=1, do_not_log=just_read) + + if (new_sim) & call ALE_regrid_accelerated(ALE_CSp, G, GV, h, tv, regrid_iterations, h, u, v) - endif endif + endif + ! This is the end of the block of code that might have initialized fields + ! internally at the start of a new run. - else ! Previous block for new_sim=.T., this block restores state -! This line calls a subroutine that reads the initial conditions ! -! from a previously generated file. ! + if (.not.new_sim) then ! This block restores the state from a restart file. + ! This line calls a subroutine that reads the initial conditions ! + ! from a previously generated file. ! call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, & G, restart_CS) if (present(Time_in)) Time = Time_in @@ -430,10 +474,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & "The exact location and properties of those sponges are \n"//& "specified via SPONGE_CONFIG.", default=.false.) if ( use_sponge ) then -! The 3 arguments here are (1) a flag indicating whether the sponge ! -! values are to be read from files, (2) the name of a file containing! -! the state toward which the model is damped, and (3) the file in ! -! which the 2-D damping rate field can be found. ! call get_param(PF, mod, "SPONGE_CONFIG", config, & "A string that sets how the sponges are configured: \n"//& " \t file - read sponge properties from the file \n"//& @@ -516,12 +556,20 @@ end subroutine MOM_initialize_state ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_thickness_from_file(h, G, GV, param_file, file_has_thickness) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Layer thicknesses, in m - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - logical, intent(in) :: file_has_thickness +!> This subroutine reads the layer thicknesses or interface heights from a file. +subroutine initialize_thickness_from_file(h, G, GV, param_file, file_has_thickness, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: file_has_thickness !< If true, this file contains layer + !! thicknesses; otherwise it contains + !! interface heights. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Arguments: h - The thickness that is being initialized. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. @@ -536,32 +584,39 @@ subroutine initialize_thickness_from_file(h, G, GV, param_file, file_has_thickne real :: dilate ! The amount by which each layer is dilated to agree ! with the bottom depth and free surface height, nondim. logical :: correct_thickness + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate character(len=40) :: mod = "initialize_thickness_from_file" ! This subroutine's name. character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") + if (.not.just_read) & + call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + + call get_param(param_file, mod, "INPUTDIR", inputdir, default=".", do_not_log=just_read) inputdir = slasher(inputdir) call get_param(param_file, mod, "THICKNESS_FILE", thickness_file, & - "The name of the thickness file.", fail_if_missing=.true.) + "The name of the thickness file.", & + fail_if_missing=.not.just_read, do_not_log=just_read) filename = trim(inputdir)//trim(thickness_file) - call log_param(param_file, mod, "INPUTDIR/THICKNESS_FILE", filename) + if (.not.just_read) call log_param(param_file, mod, "INPUTDIR/THICKNESS_FILE", filename) - if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + if ((.not.just_read) .and. (.not.file_exists(filename, G%Domain))) call MOM_error(FATAL, & " initialize_thickness_from_file: Unable to open "//trim(filename)) if (file_has_thickness) then + if (just_read) return ! All run-time parameters have been read, so return. call read_data(filename,"h",h(:,:,:),domain=G%Domain%mpp_domain) else call get_param(param_file, mod, "ADJUST_THICKNESS", correct_thickness, & "If true, all mass below the bottom removed if the \n"//& "topography is shallower than the thickness input file \n"//& - "would indicate.", default=.false.) + "would indicate.", default=.false., do_not_log=just_read) + if (just_read) return ! All run-time parameters have been read, so return. call read_data(filename,"eta",eta(:,:,:),domain=G%Domain%mpp_domain) @@ -671,11 +726,15 @@ end subroutine adjustEtaToFitBathymetry ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_thickness_uniform(h, G, GV, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Layer thicknesses, in m - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters +subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: h - The thickness that is being initialized. ! (in) G - The ocean's grid structure. @@ -689,10 +748,15 @@ subroutine initialize_thickness_uniform(h, G, GV, param_file) ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! ! positive upward, in m. ! + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! This subroutine has no run-time parameters. + call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") if (G%max_depth<=0.) call MOM_error(FATAL,"initialize_thickness_uniform: "// & @@ -731,17 +795,14 @@ subroutine initialize_thickness_search end subroutine initialize_thickness_search ! ----------------------------------------------------------------------------- -subroutine convert_thickness(h, G, GV, param_file, tv) +subroutine convert_thickness(h, G, GV, tv) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Layer thicknesses, being converted from m to H (m or kg m-2) - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables ! Arguments: h - The thickness that is being initialized. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. real, dimension(SZI_(G),SZJ_(G)) :: & p_top, p_bot real :: dz_geo(SZI_(G),SZJ_(G)) ! The change in geopotential height @@ -803,12 +864,15 @@ subroutine convert_thickness(h, G, GV, param_file, tv) end subroutine convert_thickness -subroutine depress_surface(h, G, GV, param_file, tv) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables +subroutine depress_surface(h, G, GV, param_file, tv, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), & + intent(inout) :: h !< Layer thicknesses, in H (usually m or kg m-2) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: h - The thickness that is being initialized. ! (in) G - The ocean's grid structure. ! (in) GV - The ocean's vertical grid structure. @@ -825,26 +889,33 @@ subroutine depress_surface(h, G, GV, param_file, tv) character(len=40) :: mod = "depress_surface" ! This subroutine's name. character(len=200) :: inputdir, eta_srf_file ! Strings for file/path character(len=200) :: filename, eta_srf_var ! Strings for file/path + logical :: just_read ! If true, just read parameters but set nothing. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + ! Read the surface height (or pressure) from a file. call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mod, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,& "The initial condition file for the surface height.", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, & "The initial condition variable for the surface height.",& - default="SSH") + default="SSH", do_not_log=just_read) filename = trim(inputdir)//trim(eta_srf_file) - call log_param(param_file, mod, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) + if (.not.just_read) & + call log_param(param_file, mod, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) - call read_data(filename,eta_srf_var,eta_sfc,domain=G%Domain%mpp_domain) call get_param(param_file, mod, "SURFACE_HEIGHT_IC_SCALE", scale_factor, & "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into \n"//& - "units of m", units="variable", default=1.0) + "units of m", units="variable", default=1.0, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + + call read_data(filename,eta_srf_var,eta_sfc,domain=G%Domain%mpp_domain) if (scale_factor /= 1.0) then ; do j=js,je ; do i=is,ie eta_sfc(i,j) = eta_sfc(i,j) * scale_factor @@ -885,13 +956,17 @@ end subroutine depress_surface !> Adjust the layer thicknesses by cutting away the top of each model column at the depth !! where the hydrostatic pressure matches an imposed surface pressure read from file. -subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) - type(param_file_type), intent(in) :: PF !< Parameter file structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - type(ALE_CS), pointer :: ALE_CSp !< ALE control structure - type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness (H units, m or Pa) +subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h, just_read_params) + type(param_file_type), intent(in) :: PF !< Parameter file structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(ALE_CS), pointer :: ALE_CSp !< ALE control structure + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(inout) :: h !< Layer thickness (H units, m or Pa) + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables character(len=200) :: mod = "trim_for_ice" real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface (Pa) @@ -900,30 +975,37 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path real :: scale_factor, min_thickness integer :: i, j, k + logical :: just_read ! If true, just read parameters but set nothing. logical :: use_remapping type(remapping_CS), pointer :: remap_CS => NULL() + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(PF, mod, "SURFACE_PRESSURE_FILE", p_surf_file, & "The initial condition file for the surface height.", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(PF, mod, "SURFACE_PRESSURE_VAR", p_surf_var, & "The initial condition variable for the surface height.", & - units="kg m-2", default="") + units="kg m-2", default="", do_not_log=just_read) call get_param(PF, mod, "INPUTDIR", inputdir, default=".", do_not_log=.true.) filename = trim(slasher(inputdir))//trim(p_surf_file) - call log_param(PF, mod, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) + if (.not.just_read) call log_param(PF, mod, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) - call read_data(filename, p_surf_var, p_surf, domain=G%Domain%mpp_domain) call get_param(PF, mod, "SURFACE_PRESSURE_SCALE", scale_factor, & "A scaling factor to convert SURFACE_PRESSURE_VAR from\n"//& "file SURFACE_PRESSURE_FILE into a surface pressure.", & - units="file dependent", default=1.) - if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:) + units="file dependent", default=1., do_not_log=just_read) call get_param(PF, mod, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', & - units='m', default=1.e-3) + units='m', default=1.e-3, do_not_log=just_read) call get_param(PF, mod, "TRIMMING_USES_REMAPPING", use_remapping, & 'When trimming the column, also remap T and S.', & - default=.false.) + default=.false., do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + + call read_data(filename, p_surf_var, p_surf, domain=G%Domain%mpp_domain) + if (scale_factor /= 1.) p_surf(:,:) = scale_factor * p_surf(:,:) + if (use_remapping) then allocate(remap_CS) call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true.) @@ -1035,11 +1117,14 @@ subroutine cut_off_column_top(nk, tv, Rho0, G_earth, depth, min_thickness, & end subroutine cut_off_column_top ! ----------------------------------------------------------------------------- -subroutine initialize_velocity_from_file(u, v, G, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure +subroutine initialize_velocity_from_file(u, v, G, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u !< The zonal velocity that is being initialized, in m s-1 real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v !< The meridional velocity that is being initialized, in m s-1 - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to + !! parse for modelparameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: u - The zonal velocity that is being initialized. ! (out) v - The meridional velocity that is being initialized. ! (in) G - The ocean's grid structure. @@ -1048,15 +1133,20 @@ subroutine initialize_velocity_from_file(u, v, G, param_file) ! This subroutine reads the initial velocity components from file character(len=40) :: mod = "initialize_velocity_from_file" ! This subroutine's name. character(len=200) :: filename,velocity_file,inputdir ! Strings for file/path + logical :: just_read ! If true, just read parameters but set nothing. - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") call get_param(param_file, mod, "VELOCITY_FILE", velocity_file, & "The name of the velocity initial condition file.", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) + if (just_read) return ! All run-time parameters have been read, so return. + filename = trim(inputdir)//trim(velocity_file) call log_param(param_file, mod, "INPUTDIR/VELOCITY_FILE", filename) @@ -1072,11 +1162,14 @@ end subroutine initialize_velocity_from_file ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_velocity_zero(u, v, G, param_file) +subroutine initialize_velocity_zero(u, v, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u !< The zonal velocity that is being initialized, in m s-1 real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v !< The meridional velocity that is being initialized, in m s-1 - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to + !! parse for modelparameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: u - The zonal velocity that is being initialized. ! (out) v - The meridional velocity that is being initialized. ! (in) G - The ocean's grid structure. @@ -1084,11 +1177,16 @@ subroutine initialize_velocity_zero(u, v, G, param_file) ! This subroutine sets the initial velocity components to zero character(len=200) :: mod = "initialize_velocity_zero" ! This subroutine's name. + logical :: just_read ! If true, just read parameters but set nothing. 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 ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + + if (just_read) return ! All run-time parameters have been read, so return. do k=1,nz ; do j=js,je ; do I=Isq,Ieq u(I,j,k) = 0.0 @@ -1102,11 +1200,14 @@ end subroutine initialize_velocity_zero ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_velocity_uniform(u, v, G, param_file) +subroutine initialize_velocity_uniform(u, v, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u !< The zonal velocity that is being initialized, in m s-1 real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v !< The meridional velocity that is being initialized, in m s-1 - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to + !! parse for modelparameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: u - The zonal velocity that is being initialized. ! (out) v - The meridional velocity that is being initialized. ! (in) G - The ocean's grid structure. @@ -1115,16 +1216,21 @@ subroutine initialize_velocity_uniform(u, v, G, param_file) ! This subroutine sets the initial velocity components to uniform integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz real :: initial_u_const, initial_v_const + logical :: just_read ! If true, just read parameters but set nothing. character(len=200) :: mod = "initialize_velocity_uniform" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file, mod, "INITIAL_U_CONST", initial_u_const, & "A initial uniform value for the zonal flow.", & - units="m s-1", fail_if_missing=.true.) + units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "INITIAL_V_CONST", initial_v_const, & "A initial uniform value for the meridional flow.", & - units="m s-1", fail_if_missing=.true.) + units="m s-1", fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. do k=1,nz ; do j=js,je ; do I=Isq,Ieq u(I,j,k) = initial_u_const @@ -1137,11 +1243,14 @@ end subroutine initialize_velocity_uniform ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_velocity_circular(u, v, G, param_file) +subroutine initialize_velocity_circular(u, v, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJ_(G), SZK_(G)), intent(out) :: u !< The zonal velocity that is being initialized, in m s-1 real, dimension(SZI_(G),SZJB_(G), SZK_(G)), intent(out) :: v !< The meridional velocity that is being initialized, in m s-1 - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to + !! parse for modelparameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Arguments: u - The zonal velocity that is being initialized. ! (out) v - The meridional velocity that is being initialized. ! (in) G - The ocean's grid structure. @@ -1152,14 +1261,20 @@ subroutine initialize_velocity_circular(u, v, G, param_file) character(len=200) :: mod = "initialize_velocity_circular" real :: circular_max_u real :: dpi, psi1, psi2 + logical :: just_read ! If true, just read parameters but set nothing. 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 ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file, mod, "CIRCULAR_MAX_U", circular_max_u, & "The amplitude of zonal flow from which to scale the\n"// & "circular stream function (m/s).", & - units="m s-1", default=0.) + units="m s-1", default=0., do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + dpi=acos(0.0)*2.0 ! pi do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -1190,10 +1305,12 @@ end subroutine initialize_velocity_circular ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_temp_salt_from_file(T, S, G, param_file) +subroutine initialize_temp_salt_from_file(T, S, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! This function puts the initial layer temperatures and salinities ! ! into T(:,:,:) and S(:,:,:). ! @@ -1205,49 +1322,59 @@ subroutine initialize_temp_salt_from_file(T, S, G, param_file) ! (in) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for ! model parameter values. - character(len=200) :: filename, ts_file, salt_file, inputdir ! Strings for file/path + logical :: just_read ! If true, just read parameters but set nothing. + character(len=200) :: filename, salt_filename ! Full paths to input files + character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path character(len=40) :: mod = "initialize_temp_salt_from_file" - character(len=64) :: temp_var, salt_var + character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") call get_param(param_file, mod, "TS_FILE", ts_file, & "The initial condition file for temperature.", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) filename = trim(inputdir)//trim(ts_file) - call log_param(param_file, mod, "INPUTDIR/TS_FILE", filename) - if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & - " initialize_temp_salt_from_file: Unable to open "//trim(filename)) - + if (.not.just_read) call log_param(param_file, mod, "INPUTDIR/TS_FILE", filename) call get_param(param_file, mod, "TEMP_IC_VAR", temp_var, & "The initial condition variable for potential temperature.", & - default="PTEMP") + default="PTEMP", do_not_log=just_read) call get_param(param_file, mod, "SALT_IC_VAR", salt_var, & - "The initial condition variable for salinity.", default="SALT") + "The initial condition variable for salinity.", & + default="SALT", do_not_log=just_read) + call get_param(param_file, mod, "SALT_FILE", salt_file, & + "The initial condition file for salinity.", & + default=trim(ts_file), do_not_log=just_read) -! Read the temperatures and salinities from a netcdf file. ! - call read_data(filename, temp_var, T(:,:,:), domain=G%Domain%mpp_domain) + if (just_read) return ! All run-time parameters have been read, so return. - call get_param(param_file, mod, "SALT_FILE", salt_file, & - "The initial condition file for salinity.", default=trim(ts_file)) - filename = trim(inputdir)//trim(ts_file) if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_temp_salt_from_file: Unable to open "//trim(filename)) - call read_data(filename, salt_var, S(:,:,:), domain=G%Domain%mpp_domain) +! Read the temperatures and salinities from netcdf files. ! + call read_data(filename, temp_var, T(:,:,:), domain=G%Domain%mpp_domain) + + salt_filename = trim(inputdir)//trim(salt_file) + if (.not.file_exists(salt_filename, G%Domain)) call MOM_error(FATAL, & + " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename)) + + call read_data(salt_filename, salt_var, S(:,:,:), domain=G%Domain%mpp_domain) call callTree_leave(trim(mod)//'()') end subroutine initialize_temp_salt_from_file ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_temp_salt_from_profile(T, S, G, param_file) +subroutine initialize_temp_salt_from_profile(T, S, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! This function puts the initial layer temperatures and salinities ! ! into T(:,:,:) and S(:,:,:). ! @@ -1261,14 +1388,20 @@ subroutine initialize_temp_salt_from_profile(T, S, G, param_file) ! model parameter values. real, dimension(SZK_(G)) :: T0, S0 integer :: i, j, k + logical :: just_read ! If true, just read parameters but set nothing. character(len=200) :: filename, ts_file, inputdir ! Strings for file/path character(len=40) :: mod = "initialize_temp_salt_from_profile" - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") call get_param(param_file, mod, "TS_FILE", ts_file, & "The file with the reference profiles for temperature \n"//& - "and salinity.", fail_if_missing=.true.) + "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + call get_param(param_file, mod, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) filename = trim(inputdir)//trim(ts_file) @@ -1290,13 +1423,15 @@ end subroutine initialize_temp_salt_from_profile ! ----------------------------------------------------------------------------- -subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref) +subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(EOS_type), pointer :: eqn_of_state real, intent(in) :: P_Ref + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! This function puts the initial layer temperatures and salinities ! ! into T(:,:,:) and S(:,:,:). ! @@ -1316,22 +1451,28 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref real :: drho_dS(SZK_(G)) ! Derivative of density with salinity in kg m-3 PSU-1. ! real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3. logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity. + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "initialize_temp_salt_fit" ! This subroutine's name. integer :: i, j, k, itt, nz nz = G%ke - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") call get_param(param_file, mod, "T_REF", T_Ref, & "A reference temperature used in initialization.", & - units="degC", fail_if_missing=.true.) + units="degC", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "S_REF", S_Ref, & "A reference salinity used in initialization.", units="PSU", & - default=35.0) + default=35.0, do_not_log=just_read) call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & "If true, accept the prescribed temperature and fit the \n"//& "salinity; otherwise take salinity and fit temperature.", & - default=.false.) + default=.false., do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + do k=1,nz pres(k) = P_Ref ; S0(k) = S_Ref T0(k) = T_Ref @@ -1376,10 +1517,13 @@ end subroutine initialize_temp_salt_fit ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine initialize_temp_salt_linear(T, S, G, param_file) +subroutine initialize_temp_salt_linear(T, S, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T, S type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! This subroutine initializes linear profiles for T and S according to ! reference surface layer salinity and temperature and a specified range. ! Note that the linear distribution is set up with respect to the layer @@ -1389,21 +1533,26 @@ subroutine initialize_temp_salt_linear(T, S, G, param_file) real :: S_top, T_top ! Reference salinity and temerature within surface layer real :: S_range, T_range ! Range of salinities and temperatures over the vertical real :: delta + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "initialize_temp_salt_linear" ! This subroutine's name. - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") call get_param(param_file, mod, "T_TOP", T_top, & "Initial temperature of the top surface.", & - units="degC", fail_if_missing=.true.) + units="degC", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "T_RANGE", T_range, & "Initial temperature difference (top-bottom).", & - units="degC", fail_if_missing=.true.) + units="degC", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "S_TOP", S_top, & "Initial salinity of the top surface.", & - units="PSU", fail_if_missing=.true.) + units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "S_RANGE", S_range, & "Initial salinity difference (top-bottom).", & - units="PSU", fail_if_missing=.true.) + units="PSU", fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. ! ! Prescribe salinity ! delta_S = S_range / ( G%ke - 1.0 ); @@ -1619,14 +1768,11 @@ end subroutine set_velocity_depth_min ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- -subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) - -! Determines the isopycnal interfaces and layer potential -! temperatures and salinities directly from a z-space file on a latitude- -! longitude grid. -! -! This subroutine was written by M. Harrison, with input from R. Hallberg. -! and A. Adcroft. +!> This subroutine determines the isopycnal or other coordinate interfaces and +!! layer potential temperatures and salinities directly from a z-space file on +!! a latitude-longitude grid. +subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, just_read_params) +! This subroutine was written by M. Harrison, with input from R. Hallberg & A. Adcroft. ! ! Arguments: ! (out) h - Layer thickness, in m. @@ -1637,18 +1783,22 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) ! (in) GV - The ocean's vertical grid structure. ! (in) PF - A structure indicating the open file to parse for ! model parameter values. -! (in) dirs - A structure containing several relevant directory paths. - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: h !< Layer thicknesses, in m - type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - type(param_file_type), intent(in) :: PF - type(directories), intent(in) :: dirs + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(out) :: h !< Layer thicknesses being initialized, in m + type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic + !! variables including temperature and salinity + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(param_file_type), intent(in) :: PF !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. character(len=200) :: filename ! The name of an input file containing temperature ! and salinity in z-space; also used for ice shelf area. - character(len=200) :: inputdir ! The directory where NetCDF input files are. + character(len=200) :: shelf_file ! The name of an input file used for ice shelf area. + character(len=200) :: inputdir ! The directory where NetCDF input filesare. character(len=200) :: mesg, area_varname, ice_shelf_file type(EOS_type), pointer :: eos => NULL() @@ -1673,19 +1823,18 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) real :: min_depth real :: dilate real :: missing_value_temp, missing_value_salt - logical :: new_sim logical :: correct_thickness character(len=40) :: potemp_var, salin_var character(len=8) :: laynum integer, parameter :: niter=10 ! number of iterations for t/s adjustment to layer density + logical :: just_read ! If true, just read parameters but set nothing. logical :: adjust_temperature = .true. ! fit t/s to target densities real, parameter :: missing_value = -1.e20 real, parameter :: temp_land_fill = 0.0, salt_land_fill = 35.0 logical :: reentrant_x, tripolar_n,dbg logical :: debug = .false. ! manually set this to true for verbose output - !data arrays real, dimension(:), allocatable :: z_edges_in, z_in, Rb real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z @@ -1722,12 +1871,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) PI_180=atan(1.0)/45. - call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") - call log_version(PF, mod, version, "") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - new_sim = .false. - if ((dirs%input_filename(1:1) == 'n') .and. & - (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. + if (.not.just_read) call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") + if (.not.just_read) call log_version(PF, mod, version, "") inputdir = "." ; call get_param(PF, mod, "INPUTDIR", inputdir) inputdir = slasher(inputdir) @@ -1746,35 +1893,62 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) call get_param(PF, mod, "TEMP_SALT_Z_INIT_FILE",filename, & "The name of the z-space input file used to initialize \n"//& "the layer thicknesses, temperatures and salinities.", & - default="temp_salt_z.nc") + default="temp_salt_z.nc", do_not_log=just_read) filename = trim(inputdir)//trim(filename) call get_param(PF, mod, "Z_INIT_FILE_PTEMP_VAR", potemp_var, & "The name of the potential temperature variable in \n"//& - "TEMP_SALT_Z_INIT_FILE.", default="ptemp") + "TEMP_SALT_Z_INIT_FILE.", default="ptemp", do_not_log=just_read) call get_param(PF, mod, "Z_INIT_FILE_SALT_VAR", salin_var, & "The name of the salinity variable in \n"//& - "TEMP_SALT_Z_INIT_FILE.", default="salt") + "TEMP_SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read) call get_param(PF, mod, "Z_INIT_HOMOGENIZE", homogenize, & "If True, then horizontally homogenize the interpolated \n"//& - "initial conditions.", default=.false.) + "initial conditions.", default=.false., do_not_log=just_read) call get_param(PF, mod, "Z_INIT_ALE_REMAPPING", useALEremapping, & "If True, then remap straight to model coordinate from file.",& - default=.false.) + default=.false., do_not_log=just_read) call get_param(PF, mod, "Z_INIT_REMAPPING_SCHEME", remappingScheme, & "The remapping scheme to use if using Z_INIT_ALE_REMAPPING\n"//& - "is True.", default="PPM_IH4") + "is True.", default="PPM_IH4", do_not_log=just_read) call get_param(PF, mod, "Z_INIT_REMAP_GENERAL", remap_general, & "If false, only initializes to z* coordinates.\n"//& "If true, allows initialization directly to general coordinates.",& - default=.false.) + default=.false., do_not_log=just_read) call get_param(PF, mod, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, & "If false, only reconstructs profiles for valid data points.\n"//& "If true, inserts vanished layers below the valid data.",& - default=remap_general) + default=remap_general, do_not_log=just_read) call get_param(PF, mod, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, & "If false, uses the preferred remapping algorithm for initialization.\n"//& "If true, use an older, less robust algorithm for remapping.",& - default=.true.) + default=.true., do_not_log=just_read) + call get_param(PF, mod, "ICE_SHELF", use_ice_shelf, default=.false.) + if (use_ice_shelf) then + call get_param(PF, mod, "ICE_THICKNESS_FILE", ice_shelf_file, & + "The file from which the ice bathymetry and area are read.", & + fail_if_missing=.not.just_read, do_not_log=just_read) + shelf_file = trim(inputdir)//trim(ice_shelf_file) + if (.not.just_read) call log_param(PF, mod, "INPUTDIR/THICKNESS_FILE", shelf_file) + call get_param(PF, mod, "ICE_AREA_VARNAME", area_varname, & + "The name of the area variable in ICE_THICKNESS_FILE.", & + fail_if_missing=.not.just_read, do_not_log=just_read) + endif + if (.not.useALEremapping) then + call get_param(PF, mod, "ADJUST_THICKNESS", correct_thickness, & + "If true, all mass below the bottom removed if the \n"//& + "topography is shallower than the thickness input file \n"//& + "would indicate.", default=.false., do_not_log=just_read) + + call get_param(PF, mod, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, & + "If true, all the interior layers are adjusted to \n"//& + "their target densities using mostly temperature \n"//& + "This approach can be problematic, particularly in the \n"//& + "high latitudes.", default=.true., do_not_log=just_read) + endif + if (just_read) then + call cpu_clock_end(id_clock_routine) + return ! All run-time parameters have been read, so return. + endif ! Read input grid coordinates for temperature and salinity field ! in z-coordinate dataset. The file is REQUIRED to contain the @@ -1791,7 +1965,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) ! to the North/South Pole past the limits of the input data, they are extrapolated using the average ! value at the northernmost/southernmost latitude. - call horiz_interp_and_extrap_tracer(filename, potemp_var,1.0,1, & G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, tripolar_n, homogenize) @@ -1821,20 +1994,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) call pass_var(rho_z,G%Domain) ! This is needed for building an ALE grid under ice shelves - call get_param(PF, mod, "ICE_SHELF", use_ice_shelf, default=.false.) if (use_ice_shelf) then - call get_param(PF, mod, "ICE_THICKNESS_FILE", ice_shelf_file, & - "The file from which the ice bathymetry and area are read.", & - fail_if_missing=.true.) - filename = trim(inputdir)//trim(ice_shelf_file) - call log_param(PF, mod, "INPUTDIR/THICKNESS_FILE", filename) - call get_param(PF, mod, "ICE_AREA_VARNAME", area_varname, & - "The name of the area variable in ICE_THICKNESS_FILE.", & - fail_if_missing=.true.) - if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & - "MOM_temp_salt_initialize_from_Z: Unable to open "//trim(filename)) + if (.not.file_exists(shelf_file, G%Domain)) call MOM_error(FATAL, & + "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file)) - call read_data(filename,trim(area_varname),area_shelf_h,domain=G%Domain%mpp_domain) + call read_data(shelf_file,trim(area_varname),area_shelf_h,domain=G%Domain%mpp_domain) ! initialize frac_shelf_h with zeros (open water everywhere) frac_shelf_h(:,:) = 0.0 @@ -1959,17 +2123,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, Rb, G%bathyT(is:ie,js:je), & nlevs(is:ie,js:je), nkml, nkbl, min_depth) - call get_param(PF, mod, "ADJUST_THICKNESS", correct_thickness, & - "If true, all mass below the bottom removed if the \n"//& - "topography is shallower than the thickness input file \n"//& - "would indicate.", default=.false.) - - call get_param(PF, mod, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, & - "If true, all the interior layers are adjusted to \n"//& - "their target densities using mostly temperature \n"//& - "This approach can be problematic, particularly in the \n"//& - "high latitudes.", default=.true.) - if (correct_thickness) then call adjustEtaToFitBathymetry(G, GV, zi, h) else @@ -2003,31 +2156,25 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) nlevs(is:ie,js:je)) do k=1,nz - - nPoints = 0 ; tempAvg = 0. ; saltAvg = 0. - do j=js,je - do i=is,ie - if (G%mask2dT(i,j) .ge. 1.0) then - nPoints = nPoints + 1 - tempAvg = tempAvg + tv%T(i,j,k) - saltAvg =saltAvg + tv%S(i,j,k) - endif - enddo - enddo - - ! Horizontally homogenize data to produce perfectly "flat" initial conditions - if (homogenize) then - call sum_across_PEs(nPoints) - call sum_across_PEs(tempAvg) - call sum_across_PEs(saltAvg) - if (nPoints>0) then - tempAvg = tempAvg/real(nPoints) - saltAvg = saltAvg/real(nPoints) - endif - tv%T(:,:,k) = tempAvg - tv%S(:,:,k) = saltAvg - endif - + nPoints = 0 ; tempAvg = 0. ; saltAvg = 0. + do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) >= 1.0) then + nPoints = nPoints + 1 + tempAvg = tempAvg + tv%T(i,j,k) + saltAvg =saltAvg + tv%S(i,j,k) + endif ; enddo ; enddo + + ! Horizontally homogenize data to produce perfectly "flat" initial conditions + if (homogenize) then + call sum_across_PEs(nPoints) + call sum_across_PEs(tempAvg) + call sum_across_PEs(saltAvg) + if (nPoints>0) then + tempAvg = tempAvg/real(nPoints) + saltAvg = saltAvg/real(nPoints) + endif + tv%T(:,:,k) = tempAvg + tv%S(:,:,k) = saltAvg + endif enddo endif ! useALEremapping diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 4d926f828e..35e356fdd2 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -27,7 +27,7 @@ module DOME2d_initialization public DOME2d_initialize_temperature_salinity public DOME2d_initialize_sponges -character(len=40) :: mod = "DOEM2D_initialization" !< This module's name. +character(len=40) :: mod = "DOME2D_initialization" !< This module's name. contains @@ -43,7 +43,10 @@ subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth ) integer :: i, j real :: x, bay_depth, l1, l2 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay +! This include declares and sets the variable "version". +#include "version_variable.h" + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "DOME2D_SHELF_WIDTH", dome2d_width_bay, & 'Width of shelf, as fraction of domain, in 2d DOME configuration.', & units='nondim',default=0.1) @@ -82,26 +85,35 @@ subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth ) end subroutine DOME2d_initialize_topography !> Initialize thicknesses according to coordinate mode -subroutine DOME2d_initialize_thickness ( h, G, GV, param_file ) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< Layer thicknesses - type(param_file_type), intent(in) :: param_file !< Parameter file structure +subroutine DOME2d_initialize_thickness ( h, G, GV, param_file, just_read_params ) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables - real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! - ! negative because it is positive upward. ! - real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! - ! positive upward, in m. ! + real :: e0(SZK_(GV)) ! The resting interface heights, in m, usually ! + ! negative because it is positive upward. ! + real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface ! + ! positive upward, in m. ! integer :: i, j, k, is, ie, js, je, nz real :: x real :: delta_h real :: min_thickness - character(len=40) :: verticalCoordinate real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay + logical :: just_read ! If true, just read parameters but set nothing. + character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg("MOM_initialization.F90, DOME2d_initialize_thickness: setting thickness") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) & + call MOM_mesg("MOM_initialization.F90, DOME2d_initialize_thickness: setting thickness") call get_param(param_file,mod,"MIN_THICKNESS",min_thickness, & default=1.e-3, do_not_log=.true.) @@ -114,6 +126,8 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, param_file ) call get_param(param_file, mod, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, & default=0.2, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + ! WARNING: this routine specifies the interface heights so that the last layer ! is vanished, even at maximum depth. In order to have a uniform ! layer distribution, use this line of code within the loop: @@ -202,13 +216,17 @@ end subroutine DOME2d_initialize_thickness !> Initialize temperature and salinity in the 2d DOME configuration -subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, param_file, eqn_of_state) +subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, param_file, & + eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC) real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt) real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness (m or Pa) type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: eqn_of_state !< Equation of state structure + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables integer :: i, j, k, is, ie, js, je, nz real :: x; @@ -217,11 +235,14 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, param_file, eqn_ real :: S_ref, T_ref; ! Reference salinity and temperature within surface layer real :: S_range, T_range; ! Range of salinities and temperatures over the vertical real :: xi0, xi1; + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: verticalCoordinate real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=.true.) call get_param(param_file, mod, "DOME2D_SHELF_WIDTH", dome2d_width_bay, & @@ -230,10 +251,16 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, param_file, eqn_ default=0.3, do_not_log=.true.) call get_param(param_file, mod, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, & default=0.2, do_not_log=.true.) - call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity',units='1e-3',fail_if_missing=.true.) - call get_param(param_file,mod,"T_REF",T_ref,'Refernce temperature',units='C',fail_if_missing=.true.) - call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range',units='1e-3',default=2.0) - call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range',units='1e-3',default=0.0) + call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity',units='1e-3', & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"T_REF",T_ref,'Refernce temperature',units='C', & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range', & + units='1e-3', default=2.0, do_not_log=just_read) + call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range', & + units='1e-3', default=0.0, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. T(:,:,:) = 0.0 S(:,:,:) = 0.0 diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 81d301fc08..a1a9ff3cb5 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -90,22 +90,30 @@ end subroutine DOME_initialize_topography ! ----------------------------------------------------------------------------- !> This subroutine initializes layer thicknesses for the DOME experiment -subroutine DOME_initialize_thickness(h, G, GV, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open file - !! to parse for model parameter values. - - real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! +subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + real :: e0(SZK_(GV)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! - real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface ! + real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface ! ! positive upward, in m. ! + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "DOME_initialize_thickness" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! This subroutine has no run-time parameters. + call MOM_mesg(" DOME_initialization.F90, DOME_initialize_thickness: setting thickness", 5) e0(1)=0.0 diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 885b8d97f6..e5c5ec6ea5 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -151,17 +151,19 @@ end subroutine ISOMIP_initialize_topography ! ----------------------------------------------------------------------------- !> Initialization of thicknesses -subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being - !! initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers - !! to any available thermodynamic - !! fields, including eq. of state. +subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any + !! available thermodynamic fields, including + !! the eqn. of state. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! @@ -171,23 +173,34 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) real :: x real :: delta_h, rho_range real :: min_thickness, s_sur, s_bot, t_sur, t_bot, rho_sur, rho_bot + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) & + call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") - call get_param(param_file,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3) + call get_param(param_file,mod,"MIN_THICKNESS",min_thickness, & + 'Minimum layer thickness', units='m', default=1.e-3, do_not_log=just_read) call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates - call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) - call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8) - call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) - call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur, & + 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, & + 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, & + 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,& + 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT call calculate_density(t_sur,s_sur,0.0,rho_sur,tv%eqn_of_state) @@ -223,6 +236,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) enddo ; enddo case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR ) ! Initial thicknesses for z coordinates + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie eta1D(nz+1) = -1.0*G%bathyT(i,j) do k=nz,1,-1 @@ -237,6 +251,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) enddo ; enddo case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = delta_h @@ -252,7 +267,7 @@ end subroutine ISOMIP_initialize_thickness !> Initial values for temperature and salinity subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, & - eqn_of_state) + eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC) @@ -260,14 +275,17 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness (m or Pa) type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: eqn_of_state !< Equation of state structure - ! Local variables + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables integer :: i, j, k, is, ie, js, je, nz, itt real :: x, ds, dt, rho_sur, rho_bot real :: xi0, xi1, dxi, r, S_sur, T_sur, S_bot, T_bot, S_range, T_range real :: z ! vertical position in z space character(len=40) :: verticalCoordinate, density_profile real :: rho_tmp + logical :: just_read ! If true, just read parameters but set nothing. logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity. real :: T0(SZK_(G)), S0(SZK_(G)) real :: drho_dT(SZK_(G)) ! Derivative of density with temperature in kg m-3 K-1. ! @@ -278,12 +296,18 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke pres(:) = 0.0 + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur,'Temperature at the surface (interface)', default=-1.9) - call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, 'Salinity at the surface (interface)', default=33.8) - call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, 'Temperature at the bottom (interface)', default=-1.9) - call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot,'Salinity at the bottom (interface)', default=34.55) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_T_SUR",t_sur, & + 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_S_SUR", s_sur, & + 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_T_BOT", t_bot, & + 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read) + call get_param(param_file, mod, "ISOMIP_S_BOT", s_bot, & + 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state) !write (*,*)'Density in the surface layer:', rho_sur @@ -293,6 +317,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_RHO, REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_SIGMA ) + if (just_read) return ! All run-time parameters have been read, so return. + S_range = s_sur - s_bot T_range = t_sur - t_bot !write(*,*)'S_range,T_range',S_range,T_range @@ -313,19 +339,20 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & "If true, accept the prescribed temperature and fit the \n"//& "salinity; otherwise take salinity and fit temperature.", & - default=.false.) + default=.false., do_not_log=just_read) call get_param(param_file, mod, "DRHO_DS", drho_dS1, & "Partial derivative of density with salinity.", & - units="kg m-3 PSU-1", fail_if_missing=.true.) + units="kg m-3 PSU-1", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "DRHO_DT", drho_dT1, & "Partial derivative of density with temperature.", & - units="kg m-3 K-1", fail_if_missing=.true.) + units="kg m-3 K-1", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "T_REF", T_Ref, & "A reference temperature used in initialization.", & - units="degC", fail_if_missing=.true.) + units="degC", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "S_REF", S_Ref, & "A reference salinity used in initialization.", units="PSU", & - default=35.0) + default=35.0, do_not_log=just_read) + if (just_read) return ! All run-time parameters have been read, so return. !write(*,*)'read drho_dS, drho_dT', drho_dS1, drho_dT1 @@ -335,55 +362,55 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, param_file, S_range = S_range / G%max_depth ! Convert S_range into dS/dz T_range = T_range / G%max_depth ! Convert T_range into dT/dz - do j=js,je ; do i=is,ie - xi0 = 0.0; - do k = 1,nz - !T0(k) = T_Ref; S0(k) = S_Ref - xi1 = xi0 + 0.5 * h(i,j,k); - S0(k) = S_sur + S_range * xi1; - T0(k) = T_sur + T_range * xi1; - xi0 = xi0 + h(i,j,k); - !write(*,*)'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k - enddo - - call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state) - !write(*,*)'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1) - call calculate_density(T0(1),S0(1),0.,rho_guess(1),eqn_of_state) - - if (fit_salin) then - ! A first guess of the layers' salinity. - do k=nz,1,-1 - S0(k) = max(0.0, S0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dS1) + do j=js,je ; do i=is,ie + xi0 = 0.0; + do k = 1,nz + !T0(k) = T_Ref; S0(k) = S_Ref + xi1 = xi0 + 0.5 * h(i,j,k); + S0(k) = S_sur + S_range * xi1; + T0(k) = T_sur + T_range * xi1; + xi0 = xi0 + h(i,j,k); + !write(*,*)'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k enddo - ! Refine the guesses for each layer. - do itt=1,6 - call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) - call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) - do k=1,nz - S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS1) + + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state) + !write(*,*)'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1) + call calculate_density(T0(1),S0(1),0.,rho_guess(1),eqn_of_state) + + if (fit_salin) then + ! A first guess of the layers' salinity. + do k=nz,1,-1 + S0(k) = max(0.0, S0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dS1) + enddo + ! Refine the guesses for each layer. + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz + S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS1) + enddo enddo - enddo - else - ! A first guess of the layers' temperatures. - do k=nz,1,-1 - T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT1 - enddo + else + ! A first guess of the layers' temperatures. + do k=nz,1,-1 + T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT1 + enddo - do itt=1,6 - call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) - call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) - do k=1,nz + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) - enddo - enddo - endif + enddo + enddo + endif - do k=1,nz - T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) - enddo + do k=1,nz + T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) + enddo - enddo ; enddo + enddo ; enddo case default call MOM_error(FATAL,"isomip_initialize: "// & diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index bcbc55b559..c06bc19723 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -48,14 +48,15 @@ module Phillips_initialization contains !> Initialize thickness field. -subroutine Phillips_initialize_thickness(h, G, GV, param_file) +subroutine Phillips_initialize_thickness(h, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: eta0(SZK_(G)+1) ! The 1-d nominal positions of the interfaces. real :: eta_im(SZJ_(G),SZK_(G)+1) ! A temporary array for zonal-mean eta, m. @@ -63,6 +64,7 @@ subroutine Phillips_initialize_thickness(h, G, GV, param_file) ! positive upward, in m. real :: damp_rate, jet_width, jet_height, y_2 real :: half_strat, half_depth + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "Phillips_initialize_thickness" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz @@ -71,17 +73,21 @@ subroutine Phillips_initialize_thickness(h, G, GV, param_file) eta_im(:,:) = 0.0 - call log_version(param_file, mod, version) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call log_version(param_file, mod, version) call get_param(param_file, mod, "HALF_STRAT_DEPTH", half_strat, & "The maximum depth of the ocean.", units="nondim", & - default = 0.5) + default = 0.5, do_not_log=just_read) call get_param(param_file, mod, "JET_WIDTH", jet_width, & "The width of the zonal-mean jet.", units="km", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "JET_HEIGHT", jet_height, & "The interface height scale associated with the \n"//& "zonal-mean jet.", units="m", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. half_depth = G%max_depth*half_strat eta0(1) = 0.0 ; eta0(nz+1) = -G%max_depth @@ -123,37 +129,45 @@ subroutine Phillips_initialize_thickness(h, G, GV, param_file) end subroutine Phillips_initialize_thickness !> Initialize velocity fields. -subroutine Phillips_initialize_velocity(u, v, G, GV, param_file) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] - type(param_file_type), intent(in) :: param_file !< A structure indicating - !! the open file to parse for model - !! parameter values. +subroutine Phillips_initialize_velocity(u, v, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & + intent(out) :: u !< i-component of velocity [m/s] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & + intent(out) :: v !< j-component of velocity [m/s] + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to + !! parse for modelparameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: damp_rate, jet_width, jet_height, x_2, y_2 real :: velocity_amplitude, pi integer :: i, j, k, is, ie, js, je, nz, m + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "Phillips_initialize_velocity" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - u(:,:,:) = 0.0 - v(:,:,:) = 0.0 - - pi = 4.0*atan(1.0) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - call log_version(param_file, mod, version) + if (.not.just_read) call log_version(param_file, mod, version) call get_param(param_file, mod, "VELOCITY_IC_PERTURB_AMP", velocity_amplitude, & "The magnitude of the initial velocity perturbation.", & - units="m s-1", default=0.001) + units="m s-1", default=0.001, do_not_log=just_read) call get_param(param_file, mod, "JET_WIDTH", jet_width, & "The width of the zonal-mean jet.", units="km", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "JET_HEIGHT", jet_height, & "The interface height scale associated with the \n"//& "zonal-mean jet.", units="m", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + + u(:,:,:) = 0.0 + v(:,:,:) = 0.0 + + pi = 4.0*atan(1.0) ! Use thermal wind shear to give a geostrophically balanced flow. do k=nz-1,1 ; do j=js,je ; do I=is-1,ie @@ -356,19 +370,10 @@ end subroutine Phillips_initialize_topography !! in MOM_surface_forcing.F90. * !! * !! These variables are all set in the set of subroutines (in this * -!! file) USER_initialize_bottom_depth, USER_initialize_thickness, * -!! USER_initialize_velocity, USER_initialize_temperature_salinity, * -!! USER_initialize_mixed_layer_density, USER_initialize_sponges, * -!! USER_set_coord, and USER_set_ref_profile. * -!! * -!! The names of these subroutines should be self-explanatory. They * -!! start with "USER_" to indicate that they will likely have to be * -!! modified for each simulation to set the initial conditions and * -!! boundary conditions. Most of these take two arguments: an integer * -!! argument specifying whether the fields are to be calculated * -!! internally or read from a NetCDF file; and a string giving the * -!! path to that file. If the field is initialized internally, the * -!! path is ignored. * +!! file) Phillips_initialize_thickness, Phillips_initialize_velocity, * +!! Phillips_initialize_topography and Phillips_initialize_sponges * +!! that seet up fields that are specific to the Phillips instability * +!! test case. * !! * !! Macros written all in capital letters are defined in MOM_memory.h. * !! * diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index d32ad1d6a7..f72335c458 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -23,6 +23,8 @@ module Rossby_front_2d_initialization ! Private (module-wise) parameters character(len=40) :: mod = "Rossby_front_2d_initialization" !< This module's name. +! This include declares and sets the variable "version". +#include "version_variable.h" public Rossby_front_initialize_thickness public Rossby_front_initialize_temperature_salinity @@ -36,29 +38,41 @@ module Rossby_front_2d_initialization contains !> Initialization of thicknesses in 2D Rossby front test -subroutine Rossby_front_initialize_thickness(h, G, GV, param_file ) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Thickness - type(param_file_type), intent(in) :: param_file !< Parameter file handle +subroutine Rossby_front_initialize_thickness(h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. integer :: i, j, k, is, ie, js, je, nz real :: Tz, Dml, eta, stretch, h0 real :: min_thickness, T_range, dRho_dT + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg("Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + if (.not.just_read) & + call MOM_mesg("Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness") + + if (.not.just_read) call log_version(param_file, mod, version, "") ! Read parameters needed to set thickness call get_param(param_file, mod, "MIN_THICKNESS", min_thickness, & - 'Minimum layer thickness',units='m',default=1.e-3) + 'Minimum layer thickness',units='m',default=1.e-3, do_not_log=just_read) call get_param(param_file, mod, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "T_RANGE", T_range, 'Initial temperature range', units='C', default=0.0) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + call get_param(param_file, mod, "T_RANGE", T_range, 'Initial temperature range', & + units='C', default=0.0, do_not_log=just_read) call get_param(param_file, mod, "DRHO_DT", dRho_dT, default=-0.2, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + Tz = T_range / G%max_depth select case ( coordinateMode(verticalCoordinate) ) @@ -95,28 +109,39 @@ end subroutine Rossby_front_initialize_thickness !> Initialization of temperature and salinity in the Rossby front test -subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, param_file, eqn_of_state) +subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, & + param_file, eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [deg C] real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt] real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Thickness type(param_file_type), intent(in) :: param_file !< Parameter file handle type(EOS_type), pointer :: eqn_of_state !< Equation of state structure + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. integer :: i, j, k, is, ie, js, je, nz real :: T_ref, S_ref ! Reference salinity and temerature within surface layer real :: T_range ! Range of salinities and temperatures over the vertical real :: y, zc, zi, dTdz + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: verticalCoordinate real :: PI ! 3.1415926... calculated as 4*atan(1) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) - call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity',units='1e-3',fail_if_missing=.true.) - call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature',units='C',fail_if_missing=.true.) - call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range',units='C',default=0.0) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity', units='1e-3', & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature',units='C',& + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range',& + units='C', default=0.0, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. T(:,:,:) = 0.0 S(:,:,:) = S_ref @@ -136,7 +161,7 @@ end subroutine Rossby_front_initialize_temperature_salinity !> Initialization of u and v in the Rossby front test -subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) +subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] @@ -145,21 +170,29 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: y ! Non-dimensional coordinate across channel, 0..pi real :: T_range ! Range of salinities and temperatures over the vertical real :: dUdT ! Factor to convert dT/dy into dU/dz, g*alpha/f real :: dRho_dT, zi, zc, zm, f, Ty, Dml, hAtU integer :: i, j, k, is, ie, js, je, nz + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file, mod, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) - call get_param(param_file, mod, "T_RANGE", T_range, 'Initial temperature range', units='C', default=0.0) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + call get_param(param_file, mod, "T_RANGE", T_range, 'Initial temperature range', & + units='C', default=0.0, do_not_log=just_read) call get_param(param_file, mod, "DRHO_DT", dRho_dT, default=-0.2, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + v(:,:,:) = 0.0 u(:,:,:) = 0.0 diff --git a/src/user/SCM_CVmix_tests.F90 b/src/user/SCM_CVmix_tests.F90 index 3c39776a07..03d91f4ad5 100644 --- a/src/user/SCM_CVmix_tests.F90 +++ b/src/user/SCM_CVmix_tests.F90 @@ -46,13 +46,15 @@ module SCM_CVmix_tests contains !> Initializes temperature and salinity for the SCM CVmix test example -subroutine SCM_CVmix_tests_TS_init(T, S, h, G, GV, param_file) +subroutine SCM_CVmix_tests_TS_init(T, S, h, G, GV, param_file, just_read_params) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: T !< Potential tempera\ture (degC) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: S !< Salinity (psu) real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(in) :: h !< Layer thickness (m or Pa) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV!< Vertical grid structure type(param_file_type), intent(in) :: param_file !< Input parameter structure + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Local variables real :: eta(SZK_(G)+1) ! The 1-d nominal positions of the interfaces. real :: UpperLayerTempMLD !< Upper layer Temp MLD thickness (m) @@ -64,31 +66,36 @@ subroutine SCM_CVmix_tests_TS_init(T, S, h, G, GV, param_file) real :: LowerLayerdTdz !< Temp gradient in lower layer (deg C m^{-1}) real :: LowerLayerdSdz !< Salt gradient in lower layer (PPT m^{-1}) real :: zC, DZ + logical :: just_read ! If true, just read parameters but set nothing. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - call log_version(param_file, mod, version) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call log_version(param_file, mod, version) call get_param(param_file,mod,"SCM_TEMP_MLD",UpperLayerTempMLD, & - 'Initial temp mixed layer depth', units='m',default=0.0) + 'Initial temp mixed layer depth', units='m',default=0.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_SALT_MLD",UpperLayerSaltMLD, & - 'Initial salt mixed layer depth', units='m',default=0.0) + 'Initial salt mixed layer depth', units='m',default=0.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_L1_SALT",UpperLayerSalt, & - 'Layer 2 surface salinity', units='1e-3',default=35.0) + 'Layer 2 surface salinity', units='1e-3',default=35.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_L1_TEMP",UpperLayerTemp, & - 'Layer 1 surface temperature', units='C', default=20.0) + 'Layer 1 surface temperature', units='C', default=20.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_L2_SALT",LowerLayerSalt, & - 'Layer 2 surface salinity', units='1e-3',default=35.0) + 'Layer 2 surface salinity', units='1e-3',default=35.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_L2_TEMP",LowerLayerTemp, & - 'Layer 2 surface temperature', units='C', default=20.0) + 'Layer 2 surface temperature', units='C', default=20.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_L2_DTDZ",LowerLayerdTdZ, & 'Initial temperature stratification in layer 2', & - units='C/m', default=0.00) + units='C/m', default=0.00, do_not_log=just_read) call get_param(param_file,mod,"SCM_L2_DSDZ",LowerLayerdSdZ, & 'Initial salinity stratification in layer 2', & - units='PPT/m', default=0.00) + units='PPT/m', default=0.00, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie eta(1) = 0. ! Reference to surface diff --git a/src/user/SCM_idealized_hurricane.F90 b/src/user/SCM_idealized_hurricane.F90 index d383f7150e..91ca7d9051 100644 --- a/src/user/SCM_idealized_hurricane.F90 +++ b/src/user/SCM_idealized_hurricane.F90 @@ -44,32 +44,41 @@ module SCM_idealized_hurricane contains !> Initializes temperature and salinity for the SCM idealized hurricane example -subroutine SCM_idealized_hurricane_TS_init(T, S, h, G, GV, param_file) +subroutine SCM_idealized_hurricane_TS_init(T, S, h, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC) real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (psu) real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness (m or Pa) type(param_file_type), intent(in) :: param_file !< Input parameter structure + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. ! Local variables real :: eta(SZK_(G)+1) ! The 1-d nominal positions of the interfaces. real :: S_ref, SST_ref, dTdZ, MLD real :: zC + logical :: just_read ! If true, just read parameters but set nothing. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - call log_version(param_file, mod, version) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) call log_version(param_file, mod, version) call get_param(param_file,mod,"SCM_S_REF",S_ref, & - 'Reference salinity', units='1e-3',default=35.0) + 'Reference salinity', units='1e-3',default=35.0, do_not_log=just_read) call get_param(param_file,mod,"SCM_SST_REF",SST_ref, & - 'Reference surface temperature', units='C', fail_if_missing=.true.) + 'Reference surface temperature', units='C', & + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file,mod,"SCM_DTDZ",dTdZ, & 'Initial temperature stratification below mixed layer', & - units='C/m', fail_if_missing=.true.) + units='C/m', fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file,mod,"SCM_MLD",MLD, & - 'Initial mixed layer depth', units='m', fail_if_missing=.true.) + 'Initial mixed layer depth', units='m', & + fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie eta(1) = 0. ! Reference to surface diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90 index 8211d21d74..d0d955488f 100644 --- a/src/user/adjustment_initialization.F90 +++ b/src/user/adjustment_initialization.F90 @@ -54,55 +54,68 @@ module adjustment_initialization !> Initialization of thicknesses. !! This subroutine initializes the layer thicknesses to be uniform. !------------------------------------------------------------------------------ -subroutine adjustment_initialize_thickness ( h, G, GV, param_file ) +subroutine adjustment_initialize_thickness ( h, G, GV, param_file, just_read_params) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< The thickness that is being initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model parameter values. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! ! positive upward, in m. ! - integer :: i, j, k, is, ie, js, je, nz real :: x, y, yy, delta_S_strat, dSdz, delta_S, S_ref real :: min_thickness, adjustment_width, adjustment_delta, adjustment_deltaS real :: front_wave_amp, front_wave_length, front_wave_asym real :: target_values(SZK_(G)+1) + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mod = "adjustment_initialization" ! This module's name. + integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg("initialize_thickness_uniform: setting thickness") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) & + call MOM_mesg("initialize_thickness_uniform: setting thickness") ! Parameters used by main model initialization + if (.not.just_read) call log_version(param_file, mod, version, "") call get_param(param_file,mod,"S_REF",S_ref,fail_if_missing=.true.,do_not_log=.true.) call get_param(param_file,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness', & - units='m',default=1.0e-3) + units='m',default=1.0e-3, do_not_log=just_read) ! Parameters specific to this experiment configuration call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE",verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file,mod,"ADJUSTMENT_WIDTH",adjustment_width, & "Width of frontal zone", & - units="same as x,y",fail_if_missing=.true.) + units="same as x,y", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file,mod,"DELTA_S_STRAT",delta_S_strat, & "Top-to-bottom salinity difference of stratification", & - units="1e-3",fail_if_missing=.true.) + units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file,mod,"ADJUSTMENT_DELTAS",adjustment_deltaS, & "Salinity difference across front", & - units="1e-3",fail_if_missing=.true.) + units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file,mod,"FRONT_WAVE_AMP",front_wave_amp, & "Amplitude of trans-frontal wave perturbation", & - units="same as x,y",default=0.) + units="same as x,y",default=0., do_not_log=just_read) call get_param(param_file,mod,"FRONT_WAVE_LENGTH",front_wave_length, & "Wave-length of trans-frontal wave perturbation", & - units="same as x,y",default=0.) + units="same as x,y",default=0., do_not_log=just_read) call get_param(param_file,mod,"FRONT_WAVE_ASYM",front_wave_asym, & "Amplitude of frontal asymmetric perturbation", & - default=0.) + default=0., do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. ! WARNING: this routine specifies the interface heights so that the last layer ! is vanished, even at maximum depth. In order to have a uniform @@ -195,7 +208,7 @@ end subroutine adjustment_initialize_thickness !> Initialization of temperature and salinity. !------------------------------------------------------------------------------ subroutine adjustment_initialize_temperature_salinity ( T, S, h, G, param_file, & - eqn_of_state) + eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< The temperature that is being initialized. real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< The salinity that is being initialized. @@ -203,6 +216,8 @@ subroutine adjustment_initialize_temperature_salinity ( T, S, h, G, param_file, type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model parameter values. type(EOS_type), pointer :: eqn_of_state !< Equation of state. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. integer :: i, j, k, is, ie, js, je, nz real :: x, y, yy @@ -215,26 +230,39 @@ subroutine adjustment_initialize_temperature_salinity ( T, S, h, G, param_file, real :: adjustment_width, adjustment_deltaS real :: front_wave_amp, front_wave_length, front_wave_asym real :: eta1d(SZK_(G)+1) + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + ! Parameters used by main model initialization - call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity',units='1e-3',fail_if_missing=.true.) - call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature',units='C',fail_if_missing=.true.) - call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range',units='1e-3', & - default=2.0) - call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range',units='C', & - default=0.0) + call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity', units='1e-3', & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature', units='C', & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range', units='1e-3', & + default=2.0, do_not_log=just_read) + call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range', units='C', & + default=0.0, do_not_log=just_read) ! Parameters specific to this experiment configuration BUT logged in previous s/r call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE",verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) - call get_param(param_file,mod,"ADJUSTMENT_WIDTH",adjustment_width,fail_if_missing=.true.,do_not_log=.true.) - call get_param(param_file,mod,"ADJUSTMENT_DELTAS",adjustment_deltaS,fail_if_missing=.true.,do_not_log=.true.) - call get_param(param_file,mod,"DELTA_S_STRAT",delta_S_strat,fail_if_missing=.true.,do_not_log=.true.) - call get_param(param_file,mod,"FRONT_WAVE_AMP",front_wave_amp,default=0.,do_not_log=.true.) - call get_param(param_file,mod,"FRONT_WAVE_LENGTH",front_wave_length,default=0.,do_not_log=.true.) - call get_param(param_file,mod,"FRONT_WAVE_ASYM",front_wave_asym,default=0.,do_not_log=.true.) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) + call get_param(param_file,mod,"ADJUSTMENT_WIDTH", adjustment_width, & + fail_if_missing=.not.just_read, do_not_log=.true.) + call get_param(param_file,mod,"ADJUSTMENT_DELTAS", adjustment_deltaS, & + fail_if_missing=.not.just_read, do_not_log=.true.) + call get_param(param_file,mod,"DELTA_S_STRAT", delta_S_strat, & + fail_if_missing=.not.just_read, do_not_log=.true.) + call get_param(param_file,mod,"FRONT_WAVE_AMP", front_wave_amp, default=0., & + do_not_log=.true.) + call get_param(param_file,mod,"FRONT_WAVE_LENGTH",front_wave_length, & + default=0.,do_not_log=.true.) + call get_param(param_file,mod,"FRONT_WAVE_ASYM", front_wave_asym, default=0., & + do_not_log=.true.) + + if (just_read) return ! All run-time parameters have been read, so return. T(:,:,:) = 0.0 S(:,:,:) = 0.0 diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90 index 78c78f5e81..609696f2f4 100644 --- a/src/user/baroclinic_zone_initialization.F90 +++ b/src/user/baroclinic_zone_initialization.F90 @@ -20,7 +20,8 @@ module baroclinic_zone_initialization contains !> Reads the parameters unique to this module -subroutine bcz_params(G, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, delta_T, dTdx, L_zone) +subroutine bcz_params(G, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, & + delta_T, dTdx, L_zone, just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle real, intent(out) :: S_ref !< Reference salinity (ppt) @@ -32,29 +33,47 @@ subroutine bcz_params(G, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, de real, intent(out) :: delta_T !< Temperature difference across baroclinic zone (ppt) real, intent(out) :: dTdx !< Linear temperature gradient (ppt/m) real, intent(out) :: L_zone !< Width of baroclinic zone (m) + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + logical :: just_read ! If true, just read parameters but set nothing. - call log_version(param_file, mod, version, 'Initialization of an analytic baroclninic zone') + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (.not.just_read) & + call log_version(param_file, mod, version, 'Initialization of an analytic baroclinic zone') call openParameterBlock(param_file,'BCZIC') - call get_param(param_file,mod,"S_REF",S_ref,'Reference salinity',units='ppt',default=35.) - call get_param(param_file,mod,"DSDZ",dSdz,'Salinity stratification',units='ppt/m',default=0.0) - call get_param(param_file,mod,"DELTA_S",delta_S,'Salinity difference across baroclinic zone',units='ppt',default=0.0) - call get_param(param_file,mod,"DSDX",dSdx,'Meridional salinity difference',units='ppt/'//trim(G%x_axis_units),default=0.0) - call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature',units='C',default=10.) - call get_param(param_file,mod,"DTDZ",dTdz,'Temperature stratification',units='C/m',default=0.0) - call get_param(param_file,mod,"DELTA_T",delta_T,'Temperature difference across baroclinic zone',units='C',default=0.0) - call get_param(param_file,mod,"DTDX",dTdx,'Meridional temperature difference',units='C/'//trim(G%x_axis_units),default=0.0) - call get_param(param_file,mod,"L_ZONE",L_zone,'Width of baroclinic zone',units=G%x_axis_units,default=0.5*G%len_lat) + call get_param(param_file ,mod, "S_REF", S_ref, 'Reference salinity', units='ppt', & + default=35., do_not_log=just_read) + call get_param(param_file,mod,"DSDZ",dSdz,'Salinity stratification',units='ppt/m', & + default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"DELTA_S",delta_S,'Salinity difference across baroclinic zone', & + units='ppt', default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"DSDX",dSdx,'Meridional salinity difference', & + units='ppt/'//trim(G%x_axis_units), default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"T_REF",T_ref,'Reference temperature',units='C', & + default=10., do_not_log=just_read) + call get_param(param_file,mod,"DTDZ",dTdz,'Temperature stratification',units='C/m', & + default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"DELTA_T",delta_T,'Temperature difference across baroclinic zone', & + units='C', default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"DTDX",dTdx,'Meridional temperature difference', & + units='C/'//trim(G%x_axis_units), default=0.0, do_not_log=just_read) + call get_param(param_file,mod,"L_ZONE",L_zone,'Width of baroclinic zone', & + units=G%x_axis_units, default=0.5*G%len_lat, do_not_log=just_read) call closeParameterBlock(param_file) end subroutine bcz_params !> Initialization of temperature and salinity with the baroclinic zone initial conditions -subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, param_file) +subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, param_file, & + just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [deg C] real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt] real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Thickness type(param_file_type), intent(in) :: param_file !< Parameter file handle + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. integer :: i, j, k, is, ie, js, je, nz real :: T_ref, dTdz, dTdx, delta_T ! Parameters describing temperature distribution @@ -62,9 +81,14 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, param_file) real :: L_zone ! Width of baroclinic zone real :: zc, zi, x, xd, xs, y, yd, fn real :: PI ! 3.1415926... calculated as 4*atan(1) + logical :: just_read ! If true, just read parameters but set nothing. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call bcz_params(G, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, delta_T, dTdx, L_zone) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + call bcz_params(G, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, delta_T, dTdx, L_zone, just_read_params) + + if (just_read) return ! All run-time parameters have been read, so return. T(:,:,:) = 0. S(:,:,:) = 0. diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index d00e061997..b3b52a120d 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -91,30 +91,32 @@ end subroutine benchmark_initialize_topography !! by finding the depths of interfaces in a specified latitude-dependent !! temperature profile with an exponentially decaying thermocline on top of a !! linear stratification. -subroutine benchmark_initialize_thickness(h, G, GV, param_file, eqn_of_state, P_ref) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being - !! initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open - !! file to parse for model - !! parameter values. - type(EOS_type), pointer :: eqn_of_state !< integer that selects the - !! equation of state. - real, intent(in) :: P_Ref !< The coordinate-density - !! reference pressure in Pa. - - real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! - ! negative because it is positive upward. ! - real :: e_pert(SZK_(G)+1) ! Interface height perturbations, positive ! - ! upward, in m. ! - real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface ! - ! positive upward, in m. ! +subroutine benchmark_initialize_thickness(h, G, GV, param_file, eqn_of_state, & + P_ref, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + type(EOS_type), pointer :: eqn_of_state !< integer that selects the + !! equation of state. + real, intent(in) :: P_Ref !< The coordinate-density + !! reference pressure in Pa. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + real :: e0(SZK_(GV)+1) ! The resting interface heights, in m, usually ! + ! negative because it is positive upward. ! + real :: e_pert(SZK_(GV)+1) ! Interface height perturbations, positive ! + ! upward, in m. ! + real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface ! + ! positive upward, in m. ! real :: SST ! The initial sea surface temperature, in deg C. real :: T_int ! The initial temperature of an interface, in deg C. real :: ML_depth ! The specified initial mixed layer depth, in m. real :: thermocline_scale ! The e-folding scale of the thermocline, in m. - real, dimension(SZK_(G)) :: T0, pres, S0, rho_guess, drho, drho_dT, drho_dS + real, dimension(SZK_(GV)) :: T0, pres, S0, rho_guess, drho, drho_dT, drho_dS real :: a_exp ! The fraction of the overall stratification that is exponential. real :: I_ts, I_md ! Inverse lengthscales in m-1. real :: T_frac ! A ratio of the interface temperature to the range @@ -122,11 +124,16 @@ subroutine benchmark_initialize_thickness(h, G, GV, param_file, eqn_of_state, P_ real :: err, derr_dz ! The error between the profile's temperature and the ! interface temperature for a given z and its derivative. real :: pi, z + logical :: just_read character(len=40) :: mod = "benchmark_initialize_thickness" ! This subroutine's name. integer :: i, j, k, k1, is, ie, js, je, nz, itt is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! This subroutine has no run-time parameters. + call MOM_mesg(" benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5) k1 = GV%nk_rho_varies + 1 @@ -209,7 +216,7 @@ end subroutine benchmark_initialize_thickness !> This function puts the initial layer temperatures and salinities !! into T(:,:,:) and S(:,:,:). subroutine benchmark_init_temperature_salinity(T, S, G, GV, param_file, & - eqn_of_state, P_Ref) + eqn_of_state, P_Ref, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< The potential temperature @@ -223,6 +230,8 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, param_file, & !! equation of state. real, intent(in) :: P_Ref !< The coordinate-density !! reference pressure in Pa. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: T0(SZK_(G)), S0(SZK_(G)) real :: pres(SZK_(G)) ! Reference pressure in kg m-3. ! @@ -234,11 +243,16 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, param_file, & real :: PI ! 3.1415926... calculated as 4*atan(1) real :: SST ! The initial sea surface temperature, in deg C. real :: lat + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "benchmark_init_temperature_salinity" ! This subroutine's name. integer :: i, j, k, k1, is, ie, js, je, nz, itt is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! All run-time parameters have been read, so return. + k1 = GV%nk_rho_varies + 1 do k=1,nz diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index 6190bdaad2..fcede1478e 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -38,33 +38,42 @@ module circle_obcs_initialization contains !> This subroutine initializes layer thicknesses for the circle_obcs experiment. -subroutine circle_obcs_initialize_thickness(h, G, GV, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open - !! file to parse for model parameter values. - - real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! - ! negative because it is positive upward. ! - real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! - ! positive upward, in m. ! +subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + real :: e0(SZK_(GV)+1) ! The resting interface heights, in m, usually ! + ! negative because it is positive upward. ! + real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface ! + ! positive upward, in m. ! + real :: diskrad, rad, xCenter, xRadius, lonC, latC + logical :: just_read ! This include declares and sets the variable "version". #include "version_variable.h" - character(len=40) :: mod = "circle_obcs_initialize_thickness" ! This module's name. + character(len=40) :: mod = "circle_obcs_initialization" ! This module's name. integer :: i, j, k, is, ie, js, je, nz - real :: diskrad, rad, xCenter, xRadius, lonC, latC is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg(" circle_obcs_initialization.F90, circle_obcs_initialize_thickness: setting thickness", 5) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + if (.not.just_read) & + call MOM_mesg(" circle_obcs_initialization.F90, circle_obcs_initialize_thickness: setting thickness", 5) + + if (.not.just_read) call log_version(param_file, mod, version, "") ! Parameters read by cartesian grid initialization - call log_version(param_file, mod, version, "") call get_param(param_file, mod, "DISK_RADIUS", diskrad, & "The radius of the initially elevated disk in the \n"//& "circle_obcs test case.", units=G%x_axis_units, & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. do k=1,nz e0(K) = -G%max_depth * real(k-1) / real(nz) diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index 1cb86b63d7..297ff88192 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -89,28 +89,35 @@ subroutine dense_water_initialize_topography(D, G, param_file, max_depth) end subroutine dense_water_initialize_topography !> Initialize the temperature and salinity for the dense water experiment -subroutine dense_water_initialize_TS(G, GV, param_file, eqn_of_state, T, S, h) +subroutine dense_water_initialize_TS(G, GV, param_file, eqn_of_state, T, S, h, just_read_params) type(ocean_grid_type), intent(in) :: G !< Horizontal grid control structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: eqn_of_state !< EOS structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T, S !< Output state real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: mld, S_ref, S_range, T_ref real :: zi, zmid + logical :: just_read ! If true, just read parameters but set nothing. integer :: i, j, k, nz nz = GV%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file, mod, "DENSE_WATER_MLD", mld, & "Depth of unstratified mixed layer as a fraction of the water column.", & - units="nondim", default=default_mld) + units="nondim", default=default_mld, do_not_log=just_read) call get_param(param_file, mod, "S_REF", S_ref, do_not_log=.true.) call get_param(param_file, mod, "S_RANGE", S_range, do_not_log=.true.) call get_param(param_file, mod, "T_REF", T_ref, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + ! uniform temperature everywhere T(:,:,:) = T_ref diff --git a/src/user/external_gwave_initialization.F90 b/src/user/external_gwave_initialization.F90 index 6c4c407612..7b3c90b0df 100644 --- a/src/user/external_gwave_initialization.F90 +++ b/src/user/external_gwave_initialization.F90 @@ -35,13 +35,14 @@ module external_gwave_initialization ! ----------------------------------------------------------------------------- !> This subroutine initializes layer thicknesses for the external_gwave experiment. -subroutine external_gwave_initialize_thickness(h, G, param_file) +subroutine external_gwave_initialize_thickness(h, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being - !! initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! @@ -51,20 +52,30 @@ subroutine external_gwave_initialize_thickness(h, G, param_file) ! positive upward, in m. ! real :: ssh_anomaly_height ! Vertical height of ssh anomaly real :: ssh_anomaly_width ! Lateral width of anomaly + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "external_gwave_initialize_thickness" ! This subroutine's name. +! This include declares and sets the variable "version". +#include "version_variable.h" integer :: i, j, k, is, ie, js, je, nz real :: PI, Xnondim is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg(" external_gwave_initialization.F90, external_gwave_initialize_thickness: setting thickness", 5) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + if (.not.just_read) & + call MOM_mesg(" external_gwave_initialization.F90, external_gwave_initialize_thickness: setting thickness", 5) + + if (.not.just_read) call log_version(param_file, mod, version, "") call get_param(param_file, mod, "SSH_ANOMALY_HEIGHT", ssh_anomaly_height, & "The vertical displacement of the SSH anomaly. ", units="m", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "SSH_ANOMALY_WIDTH", ssh_anomaly_width, & "The lateral width of the SSH anomaly. ", units="coordinate", & - fail_if_missing=.true.) + fail_if_missing=.not.just_read, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. + PI = 4.0*atan(1.0) do j=G%jsc,G%jec ; do i=G%isc,G%iec Xnondim = (G%geoLonT(i,j)-G%west_lon-0.5*G%len_lon) / ssh_anomaly_width diff --git a/src/user/lock_exchange_initialization.F90 b/src/user/lock_exchange_initialization.F90 index 0ec3aedfc1..1c7b3fa615 100644 --- a/src/user/lock_exchange_initialization.F90 +++ b/src/user/lock_exchange_initialization.F90 @@ -37,37 +37,49 @@ module lock_exchange_initialization !> This subroutine initializes layer thicknesses for the lock_exchange experiment. ! ----------------------------------------------------------------------------- -subroutine lock_exchange_initialize_thickness(h, G, GV, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to - !! parse for model parameter values. +subroutine lock_exchange_initialize_thickness(h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. - real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! - ! negative because it is positive upward. ! - real :: e_pert(SZK_(G)) ! Interface height perturbations, positive ! - ! upward, in m. ! - real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface ! - ! positive upward, in m. ! + real :: e0(SZK_(GV)) ! The resting interface heights, in m, usually ! + ! negative because it is positive upward. ! + real :: e_pert(SZK_(GV)) ! Interface height perturbations, positive ! + ! upward, in m. ! + real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface ! + ! positive upward, in m. ! real :: front_displacement ! Vertical displacement acrodd front real :: thermocline_thickness ! Thickness of stratified region + logical :: just_read ! If true, just read parameters but set nothing. +! This include declares and sets the variable "version". +#include "version_variable.h" character(len=40) :: mod = "lock_exchange_initialize_thickness" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg(" lock_exchange_initialization.F90, lock_exchange_initialize_thickness: setting thickness", 5) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + if (.not.just_read) & + call MOM_mesg(" lock_exchange_initialization.F90, lock_exchange_initialize_thickness: setting thickness", 5) + + if (.not.just_read) call log_version(param_file, mod, version, "") call get_param(param_file, mod, "FRONT_DISPLACEMENT", front_displacement, & "The vertical displacement of interfaces across the front. \n"//& "A value larger in magnitude that MAX_DEPTH is truncated,", & - units="m", fail_if_missing=.true.) + units="m", fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mod, "THERMOCLINE_THICKNESS", thermocline_thickness, & "The thickness of the thermocline in the lock exchange \n"//& "experiment. A value of zero creates a two layer system \n"//& "with vanished layers in between the two inflated layers.", & - default=0., units="m") + default=0., units="m", do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. do j=G%jsc,G%jec ; do i=G%isc,G%iec do k=2,nz diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 647b168277..e7d768c752 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -97,14 +97,15 @@ end subroutine seamount_initialize_topography !> Initialization of thicknesses. !! This subroutine initializes the layer thicknesses to be uniform. -subroutine seamount_initialize_thickness ( h, G, GV, param_file ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being - !! initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. +subroutine seamount_initialize_thickness ( h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: e0(SZK_(G)+1) ! The resting interface heights, in m, usually ! ! negative because it is positive upward. ! @@ -115,14 +116,20 @@ subroutine seamount_initialize_thickness ( h, G, GV, param_file ) real :: delta_h real :: min_thickness, S_surf, S_range, S_ref, S_light, S_dense character(len=20) :: verticalCoordinate + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params - call get_param(param_file,mod,"MIN_THICKNESS",min_thickness,'Minimum thickness for layer',units='m',default=1.0e-3) + if (.not.just_read) & + call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + + call get_param(param_file,mod,"MIN_THICKNESS",min_thickness, & + 'Minimum thickness for layer',& + units='m', default=1.0e-3, do_not_log=just_read) call get_param(param_file,mod,"REGRIDDING_COORDINATE_MODE",verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) ! WARNING: this routine specifies the interface heights so that the last layer ! is vanished, even at maximum depth. In order to have a uniform @@ -143,6 +150,8 @@ subroutine seamount_initialize_thickness ( h, G, GV, param_file ) call get_param(param_file, mod, "S_REF", S_ref, default=35.0, do_not_log=.true.) call get_param(param_file, mod, "TS_RANGE_S_LIGHT", S_light, default = S_Ref, do_not_log=.true.) call get_param(param_file, mod, "TS_RANGE_S_DENSE", S_dense, default = S_Ref, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + do K=1,nz+1 ! Salinity of layer k is S_light + (k-1)/(nz-1) * (S_dense - S_light) ! Salinity of interface K is S_light + (K-3/2)/(nz-1) * (S_dense - S_light) @@ -169,6 +178,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, param_file ) enddo ; enddo case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie eta1D(nz+1) = -1.0*G%bathyT(i,j) do k=nz,1,-1 @@ -180,19 +190,22 @@ subroutine seamount_initialize_thickness ( h, G, GV, param_file ) h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo - enddo ; enddo + enddo ; enddo case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie delta_h = G%bathyT(i,j) / dfloat(nz) h(i,j,:) = delta_h end do ; end do + end select end subroutine seamount_initialize_thickness !> Initial values for temperature and salinity -subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file, eqn_of_state) +subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file, & + eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC) @@ -200,27 +213,35 @@ subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness (m or Pa) type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: eqn_of_state !< Equation of state structure + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + ! Local variables integer :: i, j, k, is, ie, js, je, nz, k_light real :: xi0, xi1, dxi, r, S_surf, T_surf, S_range, T_range real :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat + logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate, density_profile is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + call get_param(param_file, mod, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & - default=DEFAULT_COORDINATE_MODE) + default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file,mod,"INITIAL_DENSITY_PROFILE", density_profile, & 'Initial profile shape. Valid values are "linear", "parabolic"\n'// & - 'and "exponential".', default='linear') + 'and "exponential".', default='linear', do_not_log=just_read) call get_param(param_file,mod,"INITIAL_SSS", S_surf, & - 'Initial surface salinity', units='1e-3', default=34.) + 'Initial surface salinity', units='1e-3', default=34., do_not_log=just_read) call get_param(param_file,mod,"INITIAL_SST", T_surf, & - 'Initial surface temperature', units='C', default=0.) + 'Initial surface temperature', units='C', default=0., do_not_log=just_read) call get_param(param_file,mod,"INITIAL_S_RANGE", S_range, & - 'Initial salinity range (bottom - surface)', units='1e-3', default=2.) + 'Initial salinity range (bottom - surface)', units='1e-3', & + default=2., do_not_log=just_read) call get_param(param_file,mod,"INITIAL_T_RANGE", T_range, & - 'Initial temperature range (bottom - surface)', units='C', default=0.) + 'Initial temperature range (bottom - surface)', units='C', & + default=0., do_not_log=just_read) select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER ) ! Initial thicknesses for layer isopycnal coordinates @@ -232,6 +253,8 @@ subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file call get_param(param_file, mod, "TS_RANGE_S_LIGHT", S_light, default = S_Ref, do_not_log=.true.) call get_param(param_file, mod, "TS_RANGE_S_DENSE", S_dense, default = S_Ref, do_not_log=.true.) call get_param(param_file, mod, "TS_RANGE_RESOLN_RATIO", res_rat, default=1.0, do_not_log=.true.) + if (just_read) return ! All run-time parameters have been read, so return. + ! Emulate the T,S used in the "ts_range" coordinate configuration code k_light = GV%nk_rho_varies + 1 do j=js,je ; do i=is,ie @@ -247,6 +270,7 @@ subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file enddo ; enddo enddo case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) ! All other coordinate use FV initialization + if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie xi0 = 0.0 do k = 1,nz diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index 7e00b24b18..275687d45b 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -82,14 +82,15 @@ end subroutine sloshing_initialize_topography !! same thickness but all interfaces (except bottom and sea surface) are !! displaced according to a half-period cosine, with maximum value on the !! left and minimum value on the right. This sets off a regular sloshing motion. -subroutine sloshing_initialize_thickness ( h, G, GV, param_file ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being - !! initialized. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. +subroutine sloshing_initialize_thickness ( h, G, GV, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized, in m. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. real :: displ(SZK_(G)+1) real :: z_unif(SZK_(G)+1) @@ -101,12 +102,17 @@ subroutine sloshing_initialize_thickness ( h, G, GV, param_file ) real :: weight_z real :: x1, y1, x2, y2 real :: t + logical :: just_read ! If true, just read parameters but set nothing. integer :: n integer :: i, j, k, is, ie, js, je, nx, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! This subroutine has no run-time parameters. + deltah = G%max_depth / nz ! Define thicknesses @@ -198,7 +204,7 @@ end subroutine sloshing_initialize_thickness !! Note that the linear distribution is set up with respect to the layer !! number, not the physical position). subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, param_file, & - eqn_of_state) + eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC). real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt). @@ -207,6 +213,8 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, param_file, & !! open file to parse for model !! parameter values. type(EOS_type), pointer :: eqn_of_state !< Equation of state structure. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. integer :: i, j, k, is, ie, js, je, nz real :: delta_S, delta_T @@ -217,17 +225,27 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, param_file, & integer :: kdelta real :: deltah real :: xi0, xi1 + logical :: just_read ! If true, just read parameters but set nothing. character(len=40) :: mod = "initialize_temp_salt_linear" ! This subroutine's ! name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - call get_param(param_file,mod,"S_REF",S_ref,'Reference value for salinity',units='1e-3',fail_if_missing=.true.) - call get_param(param_file,mod,"T_REF",T_ref,'Refernce value for temperature',units='C',fail_if_missing=.true.) + + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + call get_param(param_file,mod,"S_REF",S_ref,'Reference value for salinity', & + units='1e-3', fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file,mod,"T_REF",T_ref,'Refernce value for temperature', & + units='C', fail_if_missing=.not.just_read, do_not_log=just_read) ! The default is to assume an increase by 2 for the salinity and a uniform ! temperature - call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range.',units='1e-3',default=2.0) - call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range',units='C',default=0.0) + call get_param(param_file,mod,"S_RANGE",S_range,'Initial salinity range.', & + units='1e-3', default=2.0, do_not_log=just_read) + call get_param(param_file,mod,"T_RANGE",T_range,'Initial temperature range', & + units='C', default=0.0, do_not_log=just_read) + + if (just_read) return ! All run-time parameters have been read, so return. ! Prescribe salinity !delta_S = S_range / ( G%ke - 1.0 ) diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index cca15cd23d..4334b6fc50 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -28,7 +28,8 @@ module user_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE use MOM_io, only : write_field, slasher use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE -use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S +use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N +use MOM_open_boundary, only : OBC_DIRECTION_S use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values use MOM_variables, only : thermo_var_ptrs @@ -40,8 +41,7 @@ module user_initialization public USER_set_coord, USER_initialize_topography, USER_initialize_thickness public USER_initialize_velocity, USER_init_temperature_salinity -public USER_init_mixed_layer_density, USER_initialize_sponges -public USER_set_OBC_data, USER_set_rotation +public USER_initialize_sponges, USER_set_OBC_data, USER_set_rotation logical :: first_call = .true. @@ -89,7 +89,7 @@ subroutine USER_initialize_topography(D, G, param_file, max_depth) end subroutine USER_initialize_topography !> initialize thicknesses. -subroutine USER_initialize_thickness(h, G, param_file, T) +subroutine USER_initialize_thickness(h, G, param_file, T, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h !< The thicknesses being !! initialized. @@ -97,10 +97,19 @@ subroutine USER_initialize_thickness(h, G, param_file, T) !! open file to parse for model !! parameter values. real, intent(in), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: T !< Potential temperature. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + logical :: just_read ! If true, just read parameters but set nothing. + call MOM_error(FATAL, & "USER_initialization.F90, USER_initialize_thickness: " // & "Unmodified user routine called - you must edit the routine to use it") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! All run-time parameters have been read, so return. + h(:,:,1) = 0.0 if (first_call) call write_user_log(param_file) @@ -108,17 +117,26 @@ subroutine USER_initialize_thickness(h, G, param_file, T) end subroutine USER_initialize_thickness !> initialize velocities. -subroutine USER_initialize_velocity(u, v, G, param_file) +subroutine USER_initialize_velocity(u, v, G, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] real, dimension(SZI_(G), SZJB_(G), SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + logical :: just_read ! If true, just read parameters but set nothing. + call MOM_error(FATAL, & "USER_initialization.F90, USER_initialize_velocity: " // & "Unmodified user routine called - you must edit the routine to use it") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! All run-time parameters have been read, so return. + u(:,:,1) = 0.0 v(:,:,1) = 0.0 @@ -128,7 +146,7 @@ end subroutine USER_initialize_velocity !> This function puts the initial layer temperatures and salinities !! into T(:,:,:) and S(:,:,:). -subroutine USER_init_temperature_salinity(T, S, G, param_file, eqn_of_state) +subroutine USER_init_temperature_salinity(T, S, G, param_file, eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature (degC). real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (ppt). @@ -137,11 +155,19 @@ subroutine USER_init_temperature_salinity(T, S, G, param_file, eqn_of_state) !! parameter values. type(EOS_type), pointer :: eqn_of_state !< Integer that selects the !! equation of state. + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters without changing h. + + logical :: just_read ! If true, just read parameters but set nothing. call MOM_error(FATAL, & "USER_initialization.F90, USER_init_temperature_salinity: " // & "Unmodified user routine called - you must edit the routine to use it") + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + + if (just_read) return ! All run-time parameters have been read, so return. + T(:,:,1) = 0.0 S(:,:,1) = 0.0 @@ -149,32 +175,6 @@ subroutine USER_init_temperature_salinity(T, S, G, param_file, eqn_of_state) end subroutine USER_init_temperature_salinity -!> Set initial potential density of the mixed layer. -subroutine USER_init_mixed_layer_density(Rml, G, param_file, use_temperature, & - eqn_of_state, T, S, P_Ref) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure. - real, dimension(SZI_(G), SZJ_(G), SZK_(G)), intent(out) :: Rml !< Mixed layer potential density. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - logical, intent(in) :: use_temperature !< Whether to use potential - !! temperature. - type(EOS_type), optional, pointer :: eqn_of_state !< integer that selects the - !! equation of state. - real, dimension(SZI_(G), SZJ_(G), SZK_(G)), optional, intent(in) :: T !< Model potential temperature. - real, dimension(SZI_(G), SZJ_(G), SZK_(G)), optional, intent(in) :: S !< Model salinity. - real, optional, intent(in) :: P_Ref !< The coordinate-density - !! reference pressure in Pa. - call MOM_error(FATAL, & - "USER_initialization.F90, USER_init_mixed_layer_density: " // & - "Unmodified user routine called - you must edit the routine to use it") - - Rml(:,:,1) = 0.0 - - if (first_call) call write_user_log(param_file) - -end subroutine USER_init_mixed_layer_density - !> Set up the sponges. subroutine USER_initialize_sponges(G, use_temperature, tv, param_file, CSp, h) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure.