diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000000..37d21b8ef2 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,547 @@ +# Copyright ACCESS-NRI and contributors. See the top-level LICENSE file for details. + +cmake_minimum_required(VERSION 3.18) + +#[==============================================================================[ +# Basic project definition # +#]==============================================================================] +project(MOM6 + DESCRIPTION "Modular Ocean Model 6" + HOMEPAGE_URL https://github.com/ACCESS-NRI/MOM6 + LANGUAGES C Fortran) + +#[==============================================================================[ +# Project configuration # +#]==============================================================================] +message(STATUS "Build options") + +option(MOM6_ACCESS3 "Use ACCESS3 dependencies and install ACCESS3 libraries" OFF) +option(MOM6_OPENMP "Enable OpenMP threading" OFF) +option(MOM6_ASYMMETRIC "Use MOM asymmetric memory" OFF) +message(STATUS " - MOM6_ACCESS3 ${MOM6_ACCESS3}" ) +message(STATUS " - MOM6_OPENMP ${MOM6_OPENMP}" ) +message(STATUS " - MOM6_ASYMMETRIC ${MOM6_ASYMMETRIC}") + +#[==============================================================================[ +# Project configuration # +#]==============================================================================] + +list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake) +include(GNUInstallDirs) +include(FortranLib) +include(CMakePackageConfigHelpers) + +# Fortran compiler flags +if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fbacktrace -fconvert=big-endian -ffree-line-length-none -ffixed-line-length-none -fdefault-real-8 -fdefault-double-8") + if(${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER_EQUAL 10) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch") + endif() + set(CMAKE_Fortran_FLAGS_RELEASE "-O") + set(CMAKE_Fortran_FLAGS_DEBUG "-g -Wall -Og -ffpe-trap=zero,overflow -fcheck=bounds") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -qno-opt-dynamic-align -convert big_endian -assume byterecl -ftz -traceback -assume realloc_lhs -fp-model precise -r8") + set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -debug minimal") + set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -check uninit -check bounds -check pointers -fpe0 -check noarg_temp_created") +else() + message(WARNING "Fortran compiler with ID ${CMAKE_Fortran_COMPILER_ID} will be used with CMake default options") +endif() + +# C compiler flags +if(CMAKE_C_COMPILER_ID MATCHES "GNU") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -std=gnu99") + set(CMAKE_C_FLAGS_RELEASE "-O") + set(CMAKE_C_FLAGS_DEBUG "-g -Wall -Og -fbacktrace -ffpe-trap=invalid,zero,overflow -fcheck=bounds") +elseif(CMAKE_C_COMPILER_ID MATCHES "Intel") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -traceback -qno-opt-dynamic-align -fp-model precise -std=gnu99") + set(CMAKE_C_FLAGS_RELEASE "-O2 -debug minimal") + set(CMAKE_C_FLAGS_DEBUG "-O0 -g") +else() + message(WARNING "C compiler with ID ${CMAKE_C_COMPILER_ID} will be used with CMake default options") +endif() + +# See https://github.com/ACCESS-NRI/GFDL-generic-tracers/issues/29 +add_compile_definitions(_USE_GENERIC_TRACER) # _USE_MOM6_DIAG) + +if (MOM6_ACCESS3) + add_compile_definitions(CESMCOUPLED) +endif() + +## Fortran modules path; currently this is simply set to be the include dir +set(CMAKE_INSTALL_MODULEDIR ${CMAKE_INSTALL_INCLUDEDIR} + CACHE STRING + "Fortran module installation path (Not a cmake native variable)" +) + +#[==============================================================================[ +# External packages # +#]==============================================================================] + +find_package(fms COMPONENTS R8 REQUIRED) +find_package(GFDLGTracers REQUIRED) + +if(NOT TARGET NetCDF::NetCDF_Fortran OR NOT TARGET NetCDF::NetCDF_C) #fms provides NETCDF inteface + message(FATAL_ERROR "NetCDF interface missing from FMS package") +endif() +if (NOT NetCDF_PARALLEL) + message(FATAL_ERROR "NetCDF does not have parallel I/O support!") +endif() + +if (MOM6_ACCESS3) + find_package(Access3Share REQUIRED cdeps timing share nuopc_cap_share) + if(NOT TARGET ESMF::ESMF) + message(FATAL_ERROR "ESMF interface missing from Access3Share package") + endif() +endif() + +if(MOM6_OPENMP) + find_package(OpenMP REQUIRED COMPONENTS Fortran) +endif() + +#[==============================================================================[ +# Main definitions # +#]==============================================================================] + + +### Targets + +## MOM6 library + +# set paths for the source code and config source +set(SRC "${CMAKE_SOURCE_DIR}/src") +set(CONFIG_SRC "${CMAKE_SOURCE_DIR}/config_src") + +add_fortran_library(mom6lib mod STATIC) + +target_include_directories(mom6lib PRIVATE $) + +if(MOM6_ASYMMETRIC) + set(MOM_MEMORY_DIR $) +else() + set(MOM_MEMORY_DIR $) +endif() +target_include_directories(mom6lib PRIVATE ${MOM_MEMORY_DIR}) + +target_compile_options(mom6lib PRIVATE "$<$:${fortran_compile_flags}>") + +# link libraries +target_link_libraries(mom6lib PRIVATE FMS::fms_r8 GFDLGTracers::gtracers) + +if(MOM6_ACCESS3) + target_link_libraries(mom6lib + PRIVATE Access3::nuopc_cap_share Access3::share Access3::timing Access3::cdeps-common ESMF::ESMF + ) +endif() + +target_sources(mom6lib PRIVATE + ${SRC}/ALE/coord_adapt.F90 + ${SRC}/ALE/coord_hycom.F90 + ${SRC}/ALE/coord_rho.F90 + ${SRC}/ALE/coord_sigma.F90 + ${SRC}/ALE/coord_zlike.F90 + ${SRC}/ALE/MOM_ALE.F90 + ${SRC}/ALE/MOM_hybgen_regrid.F90 + ${SRC}/ALE/MOM_hybgen_remap.F90 + ${SRC}/ALE/MOM_hybgen_unmix.F90 + ${SRC}/ALE/MOM_regridding.F90 + ${SRC}/ALE/MOM_remapping.F90 + ${SRC}/ALE/P1M_functions.F90 + ${SRC}/ALE/P3M_functions.F90 + ${SRC}/ALE/PCM_functions.F90 + ${SRC}/ALE/PLM_functions.F90 + ${SRC}/ALE/polynomial_functions.F90 + ${SRC}/ALE/PPM_functions.F90 + ${SRC}/ALE/PQM_functions.F90 + ${SRC}/ALE/Recon1d_EMPLM_CWK.F90 + ${SRC}/ALE/Recon1d_EMPLM_WA.F90 + ${SRC}/ALE/Recon1d_EMPLM_WA_poly.F90 + ${SRC}/ALE/Recon1d_EPPM_CWK.F90 + ${SRC}/ALE/Recon1d_MPLM_CWK.F90 + ${SRC}/ALE/Recon1d_MPLM_WA.F90 + ${SRC}/ALE/Recon1d_MPLM_WA_poly.F90 + ${SRC}/ALE/Recon1d_PCM.F90 + ${SRC}/ALE/Recon1d_PLM_CW.F90 + ${SRC}/ALE/Recon1d_PLM_CWK.F90 + ${SRC}/ALE/Recon1d_PLM_hybgen.F90 + ${SRC}/ALE/Recon1d_PPM_CW.F90 + ${SRC}/ALE/Recon1d_PPM_CWK.F90 + ${SRC}/ALE/Recon1d_PPM_H4_2018.F90 + ${SRC}/ALE/Recon1d_PPM_H4_2019.F90 + ${SRC}/ALE/Recon1d_PPM_hybgen.F90 + ${SRC}/ALE/Recon1d_type.F90 + ${SRC}/ALE/regrid_consts.F90 + ${SRC}/ALE/regrid_edge_values.F90 + ${SRC}/ALE/regrid_interp.F90 + ${SRC}/ALE/regrid_solvers.F90 + + ${SRC}/core/MOM_barotropic.F90 + ${SRC}/core/MOM_boundary_update.F90 + ${SRC}/core/MOM_check_scaling.F90 + ${SRC}/core/MOM_checksum_packages.F90 + ${SRC}/core/MOM_continuity.F90 + ${SRC}/core/MOM_continuity_PPM.F90 + ${SRC}/core/MOM_CoriolisAdv.F90 + ${SRC}/core/MOM_density_integrals.F90 + ${SRC}/core/MOM_dynamics_split_RK2.F90 + ${SRC}/core/MOM_dynamics_split_RK2b.F90 + ${SRC}/core/MOM_dynamics_unsplit.F90 + ${SRC}/core/MOM_dynamics_unsplit_RK2.F90 + ${SRC}/core/MOM.F90 + ${SRC}/core/MOM_grid.F90 + ${SRC}/core/MOM_interface_heights.F90 + ${SRC}/core/MOM_isopycnal_slopes.F90 + ${SRC}/core/MOM_open_boundary.F90 + ${SRC}/core/MOM_porous_barriers.F90 + ${SRC}/core/MOM_PressureForce.F90 + ${SRC}/core/MOM_PressureForce_FV.F90 + ${SRC}/core/MOM_PressureForce_Montgomery.F90 + ${SRC}/core/MOM_stoch_eos.F90 + ${SRC}/core/MOM_transcribe_grid.F90 + ${SRC}/core/MOM_unit_tests.F90 + ${SRC}/core/MOM_variables.F90 + ${SRC}/core/MOM_verticalGrid.F90 + ${SRC}/core/MOM_forcing_type.F90 + + ${SRC}/diagnostics/MOM_debugging.F90 + ${SRC}/diagnostics/MOM_diagnose_KdWork.F90 + ${SRC}/diagnostics/MOM_diagnose_MLD.F90 + ${SRC}/diagnostics/MOM_diagnostics.F90 + ${SRC}/diagnostics/MOM_harmonic_analysis.F90 + ${SRC}/diagnostics/MOM_obsolete_diagnostics.F90 + ${SRC}/diagnostics/MOM_obsolete_params.F90 + ${SRC}/diagnostics/MOM_PointAccel.F90 + ${SRC}/diagnostics/MOM_spatial_means.F90 + ${SRC}/diagnostics/MOM_sum_output.F90 + ${SRC}/diagnostics/MOM_wave_speed.F90 + + ${SRC}/equation_of_state/MOM_EOS.F90 + ${SRC}/equation_of_state/MOM_EOS_base_type.F90 + ${SRC}/equation_of_state/MOM_EOS_Jackett06.F90 + ${SRC}/equation_of_state/MOM_EOS_linear.F90 + ${SRC}/equation_of_state/MOM_EOS_Roquet_rho.F90 + ${SRC}/equation_of_state/MOM_EOS_Roquet_SpV.F90 + ${SRC}/equation_of_state/MOM_EOS_TEOS10.F90 + ${SRC}/equation_of_state/MOM_EOS_UNESCO.F90 + ${SRC}/equation_of_state/MOM_EOS_Wright.F90 + ${SRC}/equation_of_state/MOM_EOS_Wright_full.F90 + ${SRC}/equation_of_state/MOM_EOS_Wright_red.F90 + ${SRC}/equation_of_state/MOM_temperature_convert.F90 + ${SRC}/equation_of_state/MOM_TFreeze.F90 + + ${SRC}/equation_of_state/TEOS10/gsw_chem_potential_water_t_exact.f90 + ${SRC}/equation_of_state/TEOS10/gsw_ct_freezing_exact.f90 + ${SRC}/equation_of_state/TEOS10/gsw_ct_freezing_poly.f90 + ${SRC}/equation_of_state/TEOS10/gsw_ct_from_pt.f90 + ${SRC}/equation_of_state/TEOS10/gsw_ct_from_t.f90 + ${SRC}/equation_of_state/TEOS10/gsw_entropy_part.f90 + ${SRC}/equation_of_state/TEOS10/gsw_entropy_part_zerop.f90 + ${SRC}/equation_of_state/TEOS10/gsw_gibbs.f90 + ${SRC}/equation_of_state/TEOS10/gsw_gibbs_ice.f90 + ${SRC}/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_error_functions.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_freezing_poly_coefficients.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_gibbs_ice_coefficients.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_kinds.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 + ${SRC}/equation_of_state/TEOS10/gsw_mod_toolbox.f90 + ${SRC}/equation_of_state/TEOS10/gsw_pt0_from_t.f90 + ${SRC}/equation_of_state/TEOS10/gsw_pt_from_ct.f90 + ${SRC}/equation_of_state/TEOS10/gsw_pt_from_t.f90 + ${SRC}/equation_of_state/TEOS10/gsw_rho.f90 + ${SRC}/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 + ${SRC}/equation_of_state/TEOS10/gsw_rho_second_derivatives.f90 + ${SRC}/equation_of_state/TEOS10/gsw_sp_from_sr.f90 + ${SRC}/equation_of_state/TEOS10/gsw_sr_from_sp.f90 + ${SRC}/equation_of_state/TEOS10/gsw_specvol.f90 + ${SRC}/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 + ${SRC}/equation_of_state/TEOS10/gsw_specvol_second_derivatives.f90 + ${SRC}/equation_of_state/TEOS10/gsw_t_deriv_chem_potential_water_t_exact.f90 + ${SRC}/equation_of_state/TEOS10/gsw_t_freezing_exact.f90 + ${SRC}/equation_of_state/TEOS10/gsw_t_freezing_poly.f90 + ${SRC}/equation_of_state/TEOS10/gsw_t_from_ct.f90 + + ${SRC}/framework/MOM_array_transform.F90 + ${SRC}/framework/MOM_checksums.F90 + ${SRC}/framework/MOM_coms.F90 + ${SRC}/framework/MOM_cpu_clock.F90 + ${SRC}/framework/MOM_data_override.F90 + ${SRC}/framework/MOM_diag_mediator.F90 + ${SRC}/framework/MOM_diag_remap.F90 + ${SRC}/framework/MOM_document.F90 + ${SRC}/framework/MOM_domains.F90 + ${SRC}/framework/MOM_dyn_horgrid.F90 + ${SRC}/framework/MOM_ensemble_manager.F90 + ${SRC}/framework/MOM_error_handler.F90 + ${SRC}/framework/MOM_file_parser.F90 + ${SRC}/framework/MOM_get_input.F90 + ${SRC}/framework/MOM_hor_index.F90 + ${SRC}/framework/MOM_horizontal_regridding.F90 + ${SRC}/framework/MOM_interpolate.F90 + ${SRC}/framework/MOM_intrinsic_functions.F90 + ${SRC}/framework/MOM_io.F90 + ${SRC}/framework/MOM_io_file.F90 + ${SRC}/framework/MOM_memory_macros.h + ${SRC}/framework/MOM_murmur_hash.F90 + ${SRC}/framework/MOM_netcdf.F90 + ${SRC}/framework/numerical_testing_type.F90 + ${SRC}/framework/MOM_random.F90 + ${SRC}/framework/MOM_restart.F90 + ${SRC}/framework/MOM_safe_alloc.F90 + ${SRC}/framework/MOM_string_functions.F90 + ${SRC}/framework/MOM_unique_scales.F90 + ${SRC}/framework/MOM_unit_scaling.F90 + ${SRC}/framework/MOM_write_cputime.F90 + ${SRC}/framework/posix.F90 + ${SRC}/framework/MOM_coupler_types.F90 + + ${SRC}/ice_shelf/MOM_ice_shelf_diag_mediator.F90 + ${SRC}/ice_shelf/MOM_ice_shelf_dynamics.F90 + ${SRC}/ice_shelf/MOM_ice_shelf.F90 + ${SRC}/ice_shelf/MOM_ice_shelf_initialize.F90 + ${SRC}/ice_shelf/MOM_ice_shelf_state.F90 + ${SRC}/ice_shelf/MOM_marine_ice.F90 + ${SRC}/ice_shelf/user_shelf_init.F90 + + ${SRC}/initialization/MOM_coord_initialization.F90 + ${SRC}/initialization/MOM_fixed_initialization.F90 + ${SRC}/initialization/MOM_grid_initialize.F90 + ${SRC}/initialization/MOM_shared_initialization.F90 + ${SRC}/initialization/MOM_state_initialization.F90 + ${SRC}/initialization/MOM_tracer_initialization_from_Z.F90 + + ${SRC}/ocean_data_assim/MOM_oda_driver.F90 + ${SRC}/ocean_data_assim/MOM_oda_incupd.F90 + + ${SRC}/parameterizations/CVmix/cvmix_background.F90 + ${SRC}/parameterizations/CVmix/cvmix_convection.F90 + ${SRC}/parameterizations/CVmix/cvmix_ddiff.F90 + ${SRC}/parameterizations/CVmix/cvmix_kinds_and_types.F90 + ${SRC}/parameterizations/CVmix/cvmix_kpp.F90 + ${SRC}/parameterizations/CVmix/cvmix_math.F90 + ${SRC}/parameterizations/CVmix/cvmix_put_get.F90 + ${SRC}/parameterizations/CVmix/cvmix_shear.F90 + ${SRC}/parameterizations/CVmix/cvmix_tidal.F90 + ${SRC}/parameterizations/CVmix/cvmix_utils.F90 + + ${SRC}/parameterizations/lateral/MOM_hor_visc.F90 + ${SRC}/parameterizations/lateral/MOM_interface_filter.F90 + ${SRC}/parameterizations/lateral/MOM_internal_tides.F90 + ${SRC}/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 + ${SRC}/parameterizations/lateral/MOM_load_love_numbers.F90 + ${SRC}/parameterizations/lateral/MOM_MEKE.F90 + ${SRC}/parameterizations/lateral/MOM_MEKE_types.F90 + ${SRC}/parameterizations/lateral/MOM_mixed_layer_restrat.F90 + ${SRC}/parameterizations/lateral/MOM_spherical_harmonics.F90 + ${SRC}/parameterizations/lateral/MOM_streaming_filter.F90 + ${SRC}/parameterizations/lateral/MOM_self_attr_load.F90 + ${SRC}/parameterizations/lateral/MOM_thickness_diffuse.F90 + ${SRC}/parameterizations/lateral/MOM_tidal_forcing.F90 + ${SRC}/parameterizations/lateral/MOM_wave_drag.F90 + ${SRC}/parameterizations/lateral/MOM_Zanna_Bolton.F90 + + ${SRC}/parameterizations/stochastic/MOM_stochastics.F90 + + ${SRC}/parameterizations/vertical/MOM_ALE_sponge.F90 + ${SRC}/parameterizations/vertical/MOM_bkgnd_mixing.F90 + ${SRC}/parameterizations/vertical/MOM_bulk_mixed_layer.F90 + ${SRC}/parameterizations/vertical/MOM_CVMix_conv.F90 + ${SRC}/parameterizations/vertical/MOM_CVMix_ddiff.F90 + ${SRC}/parameterizations/vertical/MOM_CVMix_KPP.F90 + ${SRC}/parameterizations/vertical/MOM_CVMix_shear.F90 + ${SRC}/parameterizations/vertical/MOM_diabatic_aux.F90 + ${SRC}/parameterizations/vertical/MOM_diabatic_driver.F90 + ${SRC}/parameterizations/vertical/MOM_diapyc_energy_req.F90 + ${SRC}/parameterizations/vertical/MOM_energetic_PBL.F90 + ${SRC}/parameterizations/vertical/MOM_entrain_diffusive.F90 + ${SRC}/parameterizations/vertical/MOM_full_convection.F90 + ${SRC}/parameterizations/vertical/MOM_geothermal.F90 + ${SRC}/parameterizations/vertical/MOM_internal_tide_input.F90 + ${SRC}/parameterizations/vertical/MOM_kappa_shear.F90 + ${SRC}/parameterizations/vertical/MOM_opacity.F90 + ${SRC}/parameterizations/vertical/MOM_regularize_layers.F90 + ${SRC}/parameterizations/vertical/MOM_set_diffusivity.F90 + ${SRC}/parameterizations/vertical/MOM_set_viscosity.F90 + ${SRC}/parameterizations/vertical/MOM_sponge.F90 + ${SRC}/parameterizations/vertical/MOM_tidal_mixing.F90 + ${SRC}/parameterizations/vertical/MOM_vert_friction.F90 + + ${SRC}/tracer/advection_test_tracer.F90 + ${SRC}/tracer/boundary_impulse_tracer.F90 + ${SRC}/tracer/DOME_tracer.F90 + ${SRC}/tracer/dyed_obc_tracer.F90 + ${SRC}/tracer/dye_example.F90 + ${SRC}/tracer/ideal_age_example.F90 + ${SRC}/tracer/ISOMIP_tracer.F90 + ${SRC}/tracer/MOM_CFC_cap.F90 + ${SRC}/tracer/MOM_generic_tracer.F90 + ${SRC}/tracer/MOM_hor_bnd_diffusion.F90 + ${SRC}/tracer/MARBL_forcing_mod.F90 + ${SRC}/tracer/MARBL_tracers.F90 + ${SRC}/tracer/MOM_neutral_diffusion.F90 + ${SRC}/tracer/MOM_OCMIP2_CFC.F90 + ${SRC}/tracer/MOM_offline_aux.F90 + ${SRC}/tracer/MOM_offline_main.F90 + ${SRC}/tracer/MOM_tracer_advect.F90 + ${SRC}/tracer/MOM_tracer_advect_schemes.F90 + ${SRC}/tracer/MOM_tracer_diabatic.F90 + ${SRC}/tracer/MOM_tracer_flow_control.F90 + ${SRC}/tracer/MOM_tracer_hor_diff.F90 + ${SRC}/tracer/MOM_tracer_registry.F90 + ${SRC}/tracer/MOM_tracer_types.F90 + ${SRC}/tracer/MOM_tracer_Z_init.F90 + ${SRC}/tracer/nw2_tracers.F90 + ${SRC}/tracer/oil_tracer.F90 + ${SRC}/tracer/pseudo_salt_tracer.F90 + ${SRC}/tracer/RGC_tracer.F90 + ${SRC}/tracer/tracer_example.F90 + + ${SRC}/user/adjustment_initialization.F90 + ${SRC}/user/baroclinic_zone_initialization.F90 + ${SRC}/user/basin_builder.F90 + ${SRC}/user/benchmark_initialization.F90 + ${SRC}/user/BFB_initialization.F90 + ${SRC}/user/BFB_surface_forcing.F90 + ${SRC}/user/circle_obcs_initialization.F90 + ${SRC}/user/dense_water_initialization.F90 + ${SRC}/user/DOME2d_initialization.F90 + ${SRC}/user/DOME_initialization.F90 + ${SRC}/user/dumbbell_initialization.F90 + ${SRC}/user/dumbbell_surface_forcing.F90 + ${SRC}/user/dyed_channel_initialization.F90 + ${SRC}/user/dyed_obcs_initialization.F90 + ${SRC}/user/external_gwave_initialization.F90 + ${SRC}/user/Idealized_Hurricane.F90 + ${SRC}/user/ISOMIP_initialization.F90 + ${SRC}/user/Kelvin_initialization.F90 + ${SRC}/user/lock_exchange_initialization.F90 + ${SRC}/user/MOM_controlled_forcing.F90 + ${SRC}/user/MOM_wave_interface.F90 + ${SRC}/user/Neverworld_initialization.F90 + ${SRC}/user/Phillips_initialization.F90 + ${SRC}/user/RGC_initialization.F90 + ${SRC}/user/Rossby_front_2d_initialization.F90 + ${SRC}/user/SCM_CVMix_tests.F90 + ${SRC}/user/seamount_initialization.F90 + ${SRC}/user/shelfwave_initialization.F90 + ${SRC}/user/sloshing_initialization.F90 + ${SRC}/user/soliton_initialization.F90 + ${SRC}/user/supercritical_initialization.F90 + ${SRC}/user/tidal_bay_initialization.F90 + ${SRC}/user/user_change_diffusivity.F90 + ${SRC}/user/user_initialization.F90 + ${SRC}/user/user_revise_forcing.F90 + + ${CONFIG_SRC}/external/database_comms/MOM_database_comms.F90 + ${CONFIG_SRC}/external/database_comms/database_client_interface.F90 + + ${CONFIG_SRC}/external/drifters/MOM_particles.F90 + ${CONFIG_SRC}/external/drifters/MOM_particles_types.F90 + + ${CONFIG_SRC}/external/ODA_hooks/kdtree.f90 + ${CONFIG_SRC}/external/ODA_hooks/ocean_da_core.F90 + ${CONFIG_SRC}/external/ODA_hooks/ocean_da_types.F90 + ${CONFIG_SRC}/external/ODA_hooks/write_ocean_obs.F90 + + ${CONFIG_SRC}/external/stochastic_physics/get_stochy_pattern.F90 + ${CONFIG_SRC}/external/stochastic_physics/stochastic_physics.F90 + + ${CONFIG_SRC}/infra/FMS2/MOM_coms_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_constants.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_couplertype_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_cpu_clock_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_data_override_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_diag_manager_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_domain_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_ensemble_manager_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_error_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_interp_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_io_infra.F90 + ${CONFIG_SRC}/infra/FMS2/MOM_time_manager.F90 + + ${CONFIG_SRC}/external/MARBL/marbl_constants_mod.F90 + ${CONFIG_SRC}/external/MARBL/marbl_interface.F90 + ${CONFIG_SRC}/external/MARBL/marbl_interface_public_types.F90 + ${CONFIG_SRC}/external/MARBL/marbl_logging.F90 +) + +if (MOM6_ACCESS3) + target_sources(mom6lib PRIVATE + ${CONFIG_SRC}/drivers/nuopc_cap/mom_cap_time.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/ocn_comp_NUOPC.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/time_utils.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/mom_cap.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/mom_cap_methods.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 + ${CONFIG_SRC}/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 + ) +else() + target_sources(mom6lib PRIVATE + ${CONFIG_SRC}/drivers/solo_driver/MOM_surface_forcing.F90 + ${CONFIG_SRC}/drivers/solo_driver/MESO_surface_forcing.F90 + ${CONFIG_SRC}/drivers/solo_driver/user_surface_forcing.F90 + ) +endif() + +#[==============================================================================[ +# Install or Export # +#]==============================================================================] + +## Library +if(MOM6_ACCESS3) + set_target_properties(mom6lib PROPERTIES + OUTPUT_NAME access-mom6lib + EXPORT_NAME mom6lib + ) + install(TARGETS mom6lib + EXPORT Mom6libTargets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} COMPONENT AccessMOM6_runtime NAMELINK_COMPONENT AccessMOM6_Development + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} COMPONENT AccessMOM6_Development + ) + # Fortran module files are a special case, as currently there is no standard + # way of handling them in CMake + target_include_directories(mom6lib PUBLIC "$") + get_target_property(mom_moddir mom6lib Fortran_MODULE_DIRECTORY) + install(FILES ${mom_moddir}/ocn_comp_nuopc.mod ${mom_moddir}/mom_cap_mod.mod + DESTINATION ${CMAKE_INSTALL_MODULEDIR}/Mom6lib + COMPONENT AccessMOM6_Development + ) + install(EXPORT Mom6libTargets + FILE Mom6libTargets.cmake + NAMESPACE Access3:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/Mom6lib + ) + + configure_package_config_file( + cmake/Mom6libConfig.cmake.in + "${CMAKE_CURRENT_BINARY_DIR}/Mom6libConfig.cmake" + INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/Mom6lib + ) + + install(FILES ${CMAKE_SOURCE_DIR}/cmake/FindNetCDF.cmake ${CMAKE_CURRENT_BINARY_DIR}/Mom6libConfig.cmake + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/Mom6lib + COMPONENT AccessMOM6_Development + ) + +else() + +# executable + add_executable(MOM6 ${CONFIG_SRC}/drivers/solo_driver/MOM_driver.F90) + target_link_libraries(MOM6 PRIVATE mom6lib) + + target_include_directories(MOM6 PRIVATE + ${SRC}/framework + ${MOM_MEMORY_DIR} + ${CONFIG_SRC}/drivers/solo_driver + ) + + set_target_properties(MOM6 PROPERTIES + LINKER_LANGUAGE Fortran + OUTPUT_NAME mom6 + ) + + install(TARGETS MOM6 + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + ) +endif() diff --git a/cmake/FindNetCDF.cmake b/cmake/FindNetCDF.cmake new file mode 100644 index 0000000000..1439ae8486 --- /dev/null +++ b/cmake/FindNetCDF.cmake @@ -0,0 +1,337 @@ +# (C) Copyright 2011- ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + +# Try to find NetCDF includes and library. +# Supports static and shared libaries and allows each component to be found in sepearte prefixes. +# +# This module defines +# +# - NetCDF_FOUND - System has NetCDF +# - NetCDF_INCLUDE_DIRS - the NetCDF include directories +# - NetCDF_VERSION - the version of NetCDF +# - NetCDF_CONFIG_EXECUTABLE - the netcdf-config executable if found +# - NetCDF_PARALLEL - Boolean True if NetCDF4 has parallel IO support via hdf5 and/or pnetcdf +# - NetCDF_HAS_PNETCDF - Boolean True if NetCDF4 has pnetcdf support +# +# Deprecated Defines +# - NetCDF_LIBRARIES - [Deprecated] Use NetCDF::NetCDF_ targets instead. +# +# +# Following components are available: +# +# - C - C interface to NetCDF (netcdf) +# - CXX - CXX4 interface to NetCDF (netcdf_c++4) +# - Fortran - Fortran interface to NetCDF (netcdff) +# +# For each component the following are defined: +# +# - NetCDF__FOUND - whether the component is found +# - NetCDF__LIBRARIES - the libraries for the component +# - NetCDF__LIBRARY_SHARED - Boolean is true if libraries for component are shared +# - NetCDF__INCLUDE_DIRS - the include directories for specified component +# - NetCDF::NetCDF_ - target of component to be used with target_link_libraries() +# +# The following paths will be searched in order if set in CMake (first priority) or environment (second priority) +# +# - NetCDF_ROOT - root of NetCDF installation +# - NetCDF_PATH - root of NetCDF installation +# +# The search process begins with locating NetCDF Include headers. If these are in a non-standard location, +# set one of the following CMake or environment variables to point to the location: +# +# - NetCDF_INCLUDE_DIR or NetCDF_${comp}_INCLUDE_DIR +# - NetCDF_INCLUDE_DIRS or NetCDF_${comp}_INCLUDE_DIR +# +# Notes: +# +# - Use "NetCDF::NetCDF_" targets only. NetCDF_LIBRARIES exists for backwards compatibility and should not be used. +# - These targets have all the knowledge of include directories and library search directories, and a single +# call to target_link_libraries will provide all these transitive properties to your target. Normally all that is +# needed to build and link against NetCDF is, e.g.: +# target_link_libraries(my_c_tgt PUBLIC NetCDF::NetCDF_C) +# - "NetCDF" is always the preferred naming for this package, its targets, variables, and environment variables +# - For compatibility, some variables are also set/checked using alternate names NetCDF4, NETCDF, or NETCDF4 +# - Environments relying on these older environment variable names should move to using a "NetCDF_ROOT" environment variable +# - Preferred component capitalization follows the CMake LANGUAGES variables: i.e., C, Fortran, CXX +# - For compatibility, alternate capitalizations are supported but should not be used. +# - If no components are defined, all components will be searched +# + +list( APPEND _possible_components C CXX Fortran ) + +## Include names for each component +set( NetCDF_C_INCLUDE_NAME netcdf.h ) +set( NetCDF_CXX_INCLUDE_NAME netcdf ) +set( NetCDF_Fortran_INCLUDE_NAME netcdf.mod ) + +## Library names for each component +set( NetCDF_C_LIBRARY_NAME netcdf ) +set( NetCDF_CXX_LIBRARY_NAME netcdf_c++4 ) +set( NetCDF_Fortran_LIBRARY_NAME netcdff ) + +## Enumerate search components +foreach( _comp ${_possible_components} ) + string( TOUPPER "${_comp}" _COMP ) + set( _arg_${_COMP} ${_comp} ) + set( _name_${_COMP} ${_comp} ) +endforeach() + +set( _search_components C) +foreach( _comp ${${CMAKE_FIND_PACKAGE_NAME}_FIND_COMPONENTS} ) + string( TOUPPER "${_comp}" _COMP ) + set( _arg_${_COMP} ${_comp} ) + list( APPEND _search_components ${_name_${_COMP}} ) + if( NOT _name_${_COMP} ) + message(SEND_ERROR "Find${CMAKE_FIND_PACKAGE_NAME}: COMPONENT ${_comp} is not a valid component. Valid components: ${_possible_components}" ) + endif() +endforeach() +list( REMOVE_DUPLICATES _search_components ) + +## Search hints for finding include directories and libraries +foreach( _comp IN ITEMS "_" "_C_" "_Fortran_" "_CXX_" ) + foreach( _name IN ITEMS NetCDF4 NetCDF NETCDF4 NETCDF ) + foreach( _var IN ITEMS ROOT PATH ) + list(APPEND _search_hints ${${_name}${_comp}${_var}} $ENV{${_name}${_comp}${_var}} ) + list(APPEND _include_search_hints + ${${_name}${_comp}INCLUDE_DIR} $ENV{${_name}${_comp}INCLUDE_DIR} + ${${_name}${_comp}INCLUDE_DIRS} $ENV{${_name}${_comp}INCLUDE_DIRS} ) + endforeach() + endforeach() +endforeach() +#Old-school HPC module env variable names +foreach( _name IN ITEMS NetCDF4 NetCDF NETCDF4 NETCDF ) + foreach( _comp IN ITEMS "_C" "_Fortran" "_CXX" ) + list(APPEND _search_hints ${${_name}} $ENV{${_name}}) + list(APPEND _search_hints ${${_name}${_comp}} $ENV{${_name}${_comp}}) + endforeach() +endforeach() + +## Find headers for each component +set(NetCDF_INCLUDE_DIRS) +set(_new_search_components) +foreach( _comp IN LISTS _search_components ) + if(NOT ${PROJECT_NAME}_NetCDF_${_comp}_FOUND) + list(APPEND _new_search_components ${_comp}) + endif() + find_file(NetCDF_${_comp}_INCLUDE_FILE + NAMES ${NetCDF_${_comp}_INCLUDE_NAME} + DOC "NetCDF ${_comp} include directory" + HINTS ${_include_search_hints} ${_search_hints} + PATH_SUFFIXES include include/netcdf + ) + mark_as_advanced(NetCDF_${_comp}_INCLUDE_FILE) + message(DEBUG "NetCDF_${_comp}_INCLUDE_FILE: ${NetCDF_${_comp}_INCLUDE_FILE}") + if( NetCDF_${_comp}_INCLUDE_FILE ) + get_filename_component(NetCDF_${_comp}_INCLUDE_FILE ${NetCDF_${_comp}_INCLUDE_FILE} ABSOLUTE) + get_filename_component(NetCDF_${_comp}_INCLUDE_DIR ${NetCDF_${_comp}_INCLUDE_FILE} DIRECTORY) + list(APPEND NetCDF_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIR}) + endif() +endforeach() +if(NetCDF_INCLUDE_DIRS) + list(REMOVE_DUPLICATES NetCDF_INCLUDE_DIRS) +endif() +set(NetCDF_INCLUDE_DIRS "${NetCDF_INCLUDE_DIRS}" CACHE STRING "NetCDF Include directory paths" FORCE) + +## Find n*-config executables for search components +foreach( _comp IN LISTS _search_components ) + if( _comp MATCHES "^(C)$" ) + set(_conf "c") + elseif( _comp MATCHES "^(Fortran)$" ) + set(_conf "f") + elseif( _comp MATCHES "^(CXX)$" ) + set(_conf "cxx4") + endif() + find_program( NetCDF_${_comp}_CONFIG_EXECUTABLE + NAMES n${_conf}-config + HINTS ${NetCDF_INCLUDE_DIRS} ${_include_search_hints} ${_search_hints} + PATH_SUFFIXES bin Bin ../bin ../../bin + DOC "NetCDF n${_conf}-config helper" ) + message(DEBUG "NetCDF_${_comp}_CONFIG_EXECUTABLE: ${NetCDF_${_comp}_CONFIG_EXECUTABLE}") +endforeach() + +set(_C_libs_flag --libs) +set(_Fortran_libs_flag --flibs) +set(_CXX_libs_flag --libs) +set(_C_includes_flag --includedir) +set(_Fortran_includes_flag --includedir) +set(_CXX_includes_flag --includedir) +function(netcdf_config exec flag output_var) + set(${output_var} False PARENT_SCOPE) + if( exec ) + execute_process( COMMAND ${exec} ${flag} RESULT_VARIABLE _ret OUTPUT_VARIABLE _val) + if( _ret EQUAL 0 ) + string( STRIP ${_val} _val ) + set( ${output_var} ${_val} PARENT_SCOPE ) + endif() + endif() +endfunction() + +## Find libraries for each component +set( NetCDF_LIBRARIES ) +foreach( _comp IN LISTS _search_components ) + string( TOUPPER "${_comp}" _COMP ) + + find_library( NetCDF_${_comp}_LIBRARY + NAMES ${NetCDF_${_comp}_LIBRARY_NAME} + DOC "NetCDF ${_comp} library" + HINTS ${NetCDF_${_comp}_INCLUDE_DIRS} ${_search_hints} + PATH_SUFFIXES lib64 lib ../lib64 ../lib ../../lib64 ../../lib ) + mark_as_advanced( NetCDF_${_comp}_LIBRARY ) + get_filename_component(NetCDF_${_comp}_LIBRARY ${NetCDF_${_comp}_LIBRARY} ABSOLUTE) + set(NetCDF_${_comp}_LIBRARY ${NetCDF_${_comp}_LIBRARY} CACHE STRING "NetCDF ${_comp} library" FORCE) + message(DEBUG "NetCDF_${_comp}_LIBRARY: ${NetCDF_${_comp}_LIBRARY}") + + if( NetCDF_${_comp}_LIBRARY ) + if( NetCDF_${_comp}_LIBRARY MATCHES ".a$" ) + set( NetCDF_${_comp}_LIBRARY_SHARED FALSE ) + set( _library_type STATIC) + else() + list( APPEND NetCDF_LIBRARIES ${NetCDF_${_comp}_LIBRARY} ) + set( NetCDF_${_comp}_LIBRARY_SHARED TRUE ) + set( _library_type SHARED) + endif() + endif() + + #Use nc-config to set per-component LIBRARIES variable if possible + netcdf_config( ${NetCDF_${_comp}_CONFIG_EXECUTABLE} ${_${_comp}_libs_flag} _val ) + if( _val ) + set( NetCDF_${_comp}_LIBRARIES ${_val} ) + if(NOT NetCDF_${_comp}_LIBRARY_SHARED AND NOT NetCDF_${_comp}_FOUND) #Static targets should use nc_config to get a proper link line with all necessary static targets. + list( APPEND NetCDF_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) + endif() + else() + set( NetCDF_${_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARY} ) + if(NOT NetCDF_${_comp}_LIBRARY_SHARED) + message(SEND_ERROR "Unable to properly find NetCDF. Found static libraries at: ${NetCDF_${_comp}_LIBRARY} but could not run nc-config: ${NetCDF_CONFIG_EXECUTABLE}") + endif() + endif() + + #Use nc-config to set per-component INCLUDE_DIRS variable if possible + netcdf_config( ${NetCDF_${_comp}_CONFIG_EXECUTABLE} ${_${_comp}_includes_flag} _val ) + if( _val ) + string( REPLACE " " ";" _val ${_val} ) + set( NetCDF_${_comp}_INCLUDE_DIRS ${_val} ) + else() + set( NetCDF_${_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIR} ) + endif() + + if( NetCDF_${_comp}_LIBRARIES AND NetCDF_${_comp}_INCLUDE_DIRS ) + set( ${CMAKE_FIND_PACKAGE_NAME}_${_arg_${_COMP}}_FOUND TRUE ) + if (NOT TARGET NetCDF::NetCDF_${_comp}) + add_library(NetCDF::NetCDF_${_comp} ${_library_type} IMPORTED) + set_target_properties(NetCDF::NetCDF_${_comp} PROPERTIES + IMPORTED_LOCATION ${NetCDF_${_comp}_LIBRARY} + INTERFACE_INCLUDE_DIRECTORIES "${NetCDF_${_comp}_INCLUDE_DIRS}" + INTERFACE_LINK_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) + endif() + endif() +endforeach() +if(NetCDF_LIBRARIES AND NetCDF_${_comp}_LIBRARY_SHARED) + list(REMOVE_DUPLICATES NetCDF_LIBRARIES) +endif() +set(NetCDF_LIBRARIES "${NetCDF_LIBRARIES}" CACHE STRING "NetCDF library targets" FORCE) + +## Find version via netcdf-config if possible +if (NetCDF_INCLUDE_DIRS) + if( NetCDF_C_CONFIG_EXECUTABLE ) + netcdf_config( ${NetCDF_C_CONFIG_EXECUTABLE} --version _vers ) + if( _vers ) + string(REGEX REPLACE ".* ((([0-9]+)\\.)+([0-9]+)).*" "\\1" NetCDF_VERSION "${_vers}" ) + endif() + else() + foreach( _dir IN LISTS NetCDF_INCLUDE_DIRS) + if( EXISTS "${_dir}/netcdf_meta.h" ) + file(STRINGS "${_dir}/netcdf_meta.h" _netcdf_version_lines + REGEX "#define[ \t]+NC_VERSION_(MAJOR|MINOR|PATCH|NOTE)") + string(REGEX REPLACE ".*NC_VERSION_MAJOR *\([0-9]*\).*" "\\1" _netcdf_version_major "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_MINOR *\([0-9]*\).*" "\\1" _netcdf_version_minor "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_PATCH *\([0-9]*\).*" "\\1" _netcdf_version_patch "${_netcdf_version_lines}") + string(REGEX REPLACE ".*NC_VERSION_NOTE *\"\([^\"]*\)\".*" "\\1" _netcdf_version_note "${_netcdf_version_lines}") + set(NetCDF_VERSION "${_netcdf_version_major}.${_netcdf_version_minor}.${_netcdf_version_patch}${_netcdf_version_note}") + unset(_netcdf_version_major) + unset(_netcdf_version_minor) + unset(_netcdf_version_patch) + unset(_netcdf_version_note) + unset(_netcdf_version_lines) + endif() + endforeach() + endif() +endif () + +## Detect additional package properties +netcdf_config(${NetCDF_C_CONFIG_EXECUTABLE} --has-parallel4 _val) +if( NOT _val MATCHES "^(yes|no)$" ) + netcdf_config(${NetCDF_C_CONFIG_EXECUTABLE} --has-parallel _val) +endif() +if( _val MATCHES "^(yes)$" ) + set(NetCDF_PARALLEL TRUE CACHE STRING "NetCDF has parallel IO capability via pnetcdf or hdf5." FORCE) +else() + set(NetCDF_PARALLEL FALSE CACHE STRING "NetCDF has no parallel IO capability." FORCE) +endif() + +## Finalize find_package +include(FindPackageHandleStandardArgs) + +if(NOT NetCDF_FOUND OR _new_search_components) + find_package_handle_standard_args( ${CMAKE_FIND_PACKAGE_NAME} + REQUIRED_VARS NetCDF_INCLUDE_DIRS NetCDF_LIBRARIES + VERSION_VAR NetCDF_VERSION + HANDLE_COMPONENTS ) +endif() + +foreach( _comp IN LISTS _search_components ) + if( NetCDF_${_comp}_FOUND ) + #Record found components to avoid duplication in NetCDF_LIBRARIES for static libraries + set(NetCDF_${_comp}_FOUND ${NetCDF_${_comp}_FOUND} CACHE BOOL "NetCDF ${_comp} Found" FORCE) + #Set a per-package, per-component found variable to communicate between multiple calls to find_package() + set(${PROJECT_NAME}_NetCDF_${_comp}_FOUND True) + endif() +endforeach() + +if( ${CMAKE_FIND_PACKAGE_NAME}_FOUND AND NOT ${CMAKE_FIND_PACKAGE_NAME}_FIND_QUIETLY AND _new_search_components) + message( STATUS "Find${CMAKE_FIND_PACKAGE_NAME} defines targets:" ) + message( STATUS " - NetCDF_VERSION [${NetCDF_VERSION}]") + message( STATUS " - NetCDF_PARALLEL [${NetCDF_PARALLEL}]") + foreach( _comp IN LISTS _new_search_components ) + string( TOUPPER "${_comp}" _COMP ) + message( STATUS " - NetCDF_${_comp}_CONFIG_EXECUTABLE [${NetCDF_${_comp}_CONFIG_EXECUTABLE}]") + if( ${CMAKE_FIND_PACKAGE_NAME}_${_arg_${_COMP}}_FOUND ) + get_filename_component(_root ${NetCDF_${_comp}_INCLUDE_DIR}/.. ABSOLUTE) + if( NetCDF_${_comp}_LIBRARY_SHARED ) + message( STATUS " - NetCDF::NetCDF_${_comp} [SHARED] [Root: ${_root}] Lib: ${NetCDF_${_comp}_LIBRARY} ") + else() + message( STATUS " - NetCDF::NetCDF_${_comp} [STATIC] [Root: ${_root}] Lib: ${NetCDF_${_comp}_LIBRARY} ") + endif() + endif() + endforeach() +endif() + +foreach( _prefix NetCDF NetCDF4 NETCDF NETCDF4 ${CMAKE_FIND_PACKAGE_NAME} ) + set( ${_prefix}_INCLUDE_DIRS ${NetCDF_INCLUDE_DIRS} ) + set( ${_prefix}_LIBRARIES ${NetCDF_LIBRARIES}) + set( ${_prefix}_VERSION ${NetCDF_VERSION} ) + set( ${_prefix}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_FOUND} ) + set( ${_prefix}_CONFIG_EXECUTABLE ${NetCDF_CONFIG_EXECUTABLE} ) + set( ${_prefix}_PARALLEL ${NetCDF_PARALLEL} ) + + foreach( _comp ${_search_components} ) + string( TOUPPER "${_comp}" _COMP ) + set( _arg_comp ${_arg_${_COMP}} ) + set( ${_prefix}_${_comp}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) + set( ${_prefix}_${_COMP}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) + set( ${_prefix}_${_arg_comp}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) + + set( ${_prefix}_${_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) + set( ${_prefix}_${_COMP}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) + set( ${_prefix}_${_arg_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) + + set( ${_prefix}_${_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) + set( ${_prefix}_${_COMP}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) + set( ${_prefix}_${_arg_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) + endforeach() +endforeach() diff --git a/cmake/FortranLib.cmake b/cmake/FortranLib.cmake new file mode 100644 index 0000000000..6eaec978be --- /dev/null +++ b/cmake/FortranLib.cmake @@ -0,0 +1,10 @@ +# Copyright ACCESS-NRI and contributors. See the top-level LICENSE file for details. + +function(add_fortran_library LIB MOD_DIR) + add_library(${LIB} ${ARGN}) + + get_target_property(LIB_DIR ${LIB} BINARY_DIR) + set_target_properties(${LIB} PROPERTIES Fortran_MODULE_DIRECTORY ${LIB_DIR}/${MOD_DIR}) + + target_include_directories(${LIB} INTERFACE "$") + endfunction(add_fortran_library) diff --git a/cmake/Mom6libConfig.cmake.in b/cmake/Mom6libConfig.cmake.in new file mode 100644 index 0000000000..d9ee377e37 --- /dev/null +++ b/cmake/Mom6libConfig.cmake.in @@ -0,0 +1,31 @@ +# Copyright ACCESS-NRI and contributors. See the top-level LICENSE file for details. + +@PACKAGE_INIT@ + +if(NOT Mom6lib_FIND_QUIETLY) + message(STATUS "Found Mom6lib: ${PACKAGE_PREFIX_DIR}") +endif() + +include(CMakeFindDependencyMacro) + +# Request components +set(_required_components ${Mom6lib_FIND_COMPONENTS}) + +find_dependency(fms COMPONENTS R8 REQUIRED) +find_dependency(GFDLGTracers REQUIRED) +find_dependency(NetCDF REQUIRED Fortran) +if (NOT NetCDF_PARALLEL) + message(FATAL_ERROR "NetCDF does not have parallel I/O support!") +endif() + +if(@OPENMP@) + find_package(OpenMP REQUIRED) +endif() + +# Run the normal Targets.cmake +list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}) +include("${CMAKE_CURRENT_LIST_DIR}/Mom6libTargets.cmake") +list(REMOVE_ITEM CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}) + +# Check the requested components are valid +check_required_components(_required_components) diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index e3b7b0cec7..4285fcda19 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -372,7 +372,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas !allocate(OS%sfc_state) call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., & - gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot) + gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot, & + use_iceshelves=OS%use_ice_shelf) if (present(wind_stagger)) then call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & diff --git a/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 index d1c46f4254..24e547b0e7 100644 --- a/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 +++ b/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 @@ -362,7 +362,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! Consider using a run-time flag to determine whether to do the diagnostic ! vertical integrals, since the related 3-d sums are not negligible in cost. call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, & - do_integrals=.true., gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot) + do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & + use_meltpot=use_melt_pot, use_iceshelves=OS%use_ice_shelf) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index e0e618c76e..13df7c82e0 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -2,8 +2,9 @@ module MOM_cap_mod +use field_manager_mod, only: field_manager_init, field_manager_end use MOM_domains, only: get_domain_extent -use MOM_io, only: stdout, io_infra_end +use MOM_io, only: stdout, io_infra_end, slasher use mpp_domains_mod, only: mpp_get_compute_domains use mpp_domains_mod, only: mpp_get_ntile_count, mpp_get_pelist, mpp_get_global_domain use mpp_domains_mod, only: mpp_get_domain_npes @@ -17,7 +18,7 @@ module MOM_cap_mod use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file use MOM_get_input, only: get_MOM_input, directories use MOM_domains, only: pass_var, pe_here -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type @@ -30,6 +31,7 @@ module MOM_cap_mod use MOM_cap_methods, only: ChkErr use MOM_ensemble_manager, only: ensemble_manager_init use MOM_coms, only: sum_across_PEs +use MOM_coupler_types, only: coupler_1d_bc_type, coupler_2d_bc_type #ifdef CESMCOUPLED use shr_log_mod, only: shr_log_setLogUnit @@ -37,6 +39,15 @@ module MOM_cap_mod #endif use time_utils_mod, only: esmf2fms_time +#ifdef _USE_GENERIC_TRACER +use MOM_coupler_types, only: coupler_type_spawn, coupler_type_destructor +use MOM_coupler_types, only: coupler_type_set_diags, coupler_type_send_data, coupler_type_data_override +use MOM_data_override, only: data_override_init, data_override +use MOM_cap_gtracer_flux, only: gas_exchange_init, gas_fields_restore, gas_fields_restart +use MOM_cap_gtracer_flux, only: get_coupled_field_name, add_gas_fluxes_param, UNKNOWN_CMEPS_FIELD +use MOM_cap_gtracer_flux, only: atmos_ocean_fluxes_calc +#endif + use, intrinsic :: iso_fortran_env, only: output_unit use ESMF, only: ESMF_ClockAdvance, ESMF_ClockGet, ESMF_ClockPrint, ESMF_VMget @@ -100,11 +111,14 @@ module MOM_cap_mod public SetServices public SetVM -!> Internal state type with pointers to three types defined by MOM. +!> Internal state type with pointers to types defined by MOM. type ocean_internalstate_type type(ocean_public_type), pointer :: ocean_public_type_ptr type(ocean_state_type), pointer :: ocean_state_type_ptr type(ice_ocean_boundary_type), pointer :: ice_ocean_boundary_type_ptr +#ifdef _USE_GENERIC_TRACER + type(coupler_2d_bc_type), pointer :: coupler_2d_bc_type_ptr +#endif end type !> Wrapper-derived type required to associate an internal state instance @@ -420,6 +434,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) type (ocean_public_type), pointer :: ocean_public => NULL() type (ocean_state_type), pointer :: ocean_state => NULL() type(ice_ocean_boundary_type), pointer :: Ice_ocean_boundary => NULL() + type(coupler_1d_bc_type), pointer :: gas_fields_atm => NULL() + type(coupler_1d_bc_type), pointer :: gas_fields_ocn => NULL() + type(coupler_1d_bc_type), pointer :: gas_fluxes => NULL() + type(coupler_2d_bc_type), pointer :: atm_fields => NULL() type(ocean_internalstate_wrapper) :: ocean_internalstate type(ocean_grid_type), pointer :: ocean_grid => NULL() type(directories) :: dirs @@ -452,6 +470,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) character(len=512) :: restartfile ! Path/Name of restart file character(len=2048) :: restartfiles ! Path/Name of restart files ! (same as restartfile if single restart file) + character(240) :: additional_restart_dir character(len=*), parameter :: subname='(MOM_cap:InitializeAdvertise)' character(len=32) :: calendar character(len=17) :: timestamp @@ -548,6 +567,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call MOM_infra_init(mpi_comm_mom) + call field_manager_init + ! determine the calendar if (cesm_coupled) then call NUOPC_CompAttributeGet(gcomp, name="calendar", value=cvalue, & @@ -706,13 +727,44 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif + ! Set NUOPC attribute additional_restart_dir to RESTART/ if not defined + additional_restart_dir = "RESTART/" + call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + additional_restart_dir = slasher(cvalue) + else + call ESMF_LogWrite('MOM_cap:additional_restart_dir unset. Defaulting to '//trim(additional_restart_dir), & + ESMF_LOGMSG_INFO) + endif + call NUOPC_CompAttributeSet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + ocean_public%is_ocean_pe = .true. +#ifdef _USE_GENERIC_TRACER + ! Initialise structures for extra tracer fluxes + call gas_exchange_init(gas_fields_atm=gas_fields_atm, gas_fields_ocn=gas_fields_ocn, gas_fluxes=gas_fluxes) + + if (cesm_coupled .and. len_trim(inst_suffix)>0) then + call ocean_model_init(ocean_public, ocean_state, time0, time_start, gas_fields_ocn=gas_fields_ocn, & + input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index) + else + call ocean_model_init(ocean_public, ocean_state, time0, time_start, gas_fields_ocn=gas_fields_ocn, & + input_restart_file=trim(adjustl(restartfiles))) + endif + + ! Enable data override via the data_table using the component name 'OCN' + call get_ocean_grid(ocean_state, ocean_grid) + call data_override_init(ocean_grid%Domain) +#else if (cesm_coupled .and. len_trim(inst_suffix)>0) then call ocean_model_init(ocean_public, ocean_state, time0, time_start, & input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index) else call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(adjustl(restartfiles))) endif +#endif ! GMM, this call is not needed in CESM. Check with EMC if it can be deleted. call ocean_model_flux_init(ocean_state) @@ -780,6 +832,31 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif endif +#ifdef _USE_GENERIC_TRACER + ! Allocate fields for extra tracer fluxes in Ice_ocean_boundary + ! Annoyingly, spawning doesn't copy param array, so add manually + call coupler_type_spawn(gas_fluxes, Ice_ocean_boundary%fluxes, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix='_ice_ocn') + call add_gas_fluxes_param(Ice_ocean_boundary%fluxes) + + ! Initialise structure for atmos fields related to extra tracer fluxes + ! This is set in the ESMF Internal State to be accessed elsewhere + ! TODO: should we deallocate atm_fields in a finalise step? Ice_ocean_boundary is handled + ! in a similar way and does not appear to be deallocated. + allocate(atm_fields) + ocean_internalstate%ptr%coupler_2d_bc_type_ptr => atm_fields + call coupler_type_spawn(gas_fields_atm, atm_fields, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix='_atm') + + ! Register diagnosics for extra tracer flux structures + call coupler_type_set_diags(Ice_ocean_boundary%fluxes, "ocean_flux", ocean_public%axes(1:2), time_start) + call coupler_type_set_diags(atm_fields, "atmos_sfc", ocean_public%axes(1:2), time_start) + + ! Restore ocean fields related to extra tracer fluxes from restart files + call get_MOM_input(dirs=dirs) + call gas_fields_restore(ocean_public%fields, ocean_public%domain, dirs%restart_input_dir) +#endif + if (use_waves) then if (wave_method == "EFACTOR") then allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec), source=0.0) @@ -881,6 +958,15 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) endif endif +#ifdef _USE_GENERIC_TRACER + ! Add import fields required for extra tracer fluxes + do n = 1, gas_fluxes%num_bcs + stdname = get_coupled_field_name(gas_fluxes%bc(n)%name) + if (stdname /= UNKNOWN_CMEPS_FIELD) & + call fld_list_add(fldsToOcn_num, fldsToOcn, stdname, "will provide") + enddo +#endif + !--------- export fields ------------- call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_omask" , "will provide") call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_t" , "will provide") @@ -895,6 +981,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Faoo_fco2_ocn", "will provide") endif + ! TODO: dts: How to handle export fields from generic tracers? + do n = 1,fldsToOcn_num call NUOPC_Advertise(importState, standardName=fldsToOcn(n)%stdname, name=fldsToOcn(n)%shortname, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1250,7 +1338,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) frmt = "('ERROR: ESMF mesh and MOM6 domain masks are inconsistent! - "//& "MOM n, maskMesh(n), mask(n) = ',3(i8,2x))" write(err_msg, frmt)n,maskMesh(n),mask(n) - call MOM_error(FATAL, err_msg) + call MOM_error(WARNING, err_msg) endif end do @@ -1706,11 +1794,14 @@ subroutine ModelAdvance(gcomp, rc) type (ocean_public_type), pointer :: ocean_public => NULL() type (ocean_state_type), pointer :: ocean_state => NULL() type(ice_ocean_boundary_type), pointer :: Ice_ocean_boundary => NULL() + type(coupler_2d_bc_type), pointer :: atm_fields => NULL() type(ocean_internalstate_wrapper) :: ocean_internalstate type(ocean_grid_type) , pointer :: ocean_grid type(time_type) :: Time + type(time_type) :: Time_import type(time_type) :: Time_step_coupled type(time_type) :: Time_restart_current + integer :: isc,iec,jsc,jec integer :: dth, dtm, dts integer :: nc type(ESMF_Time) :: MyTime @@ -1722,13 +1813,14 @@ subroutine ModelAdvance(gcomp, rc) integer :: writeunit integer :: localPet type(ESMF_VM) :: vm - integer :: n, i + integer :: m, n, i character(240) :: import_timestr, export_timestr character(len=128) :: fldname character(len=*),parameter :: subname='(MOM_cap:ModelAdvance)' character(len=8) :: suffix character(len=:), allocatable :: rpointer_filename character(len=17) :: timestamp + character(240) :: additional_restart_dir integer :: num_rest_files real(8) :: MPI_Wtime, timers logical :: write_restart, write_restartfh @@ -1769,6 +1861,7 @@ subroutine ModelAdvance(gcomp, rc) Time_step_coupled = esmf2fms_time(timeStep) Time = esmf2fms_time(currTime) + Time_import = Time !--------------- ! Apply ocean lag for startup runs: @@ -1844,8 +1937,39 @@ subroutine ModelAdvance(gcomp, rc) ! Import data !--------------- +#ifdef _USE_GENERIC_TRACER + atm_fields => ocean_internalstate%ptr%coupler_2d_bc_type_ptr + + call mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, atm_fields=atm_fields, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! Potentially override atm_fields from data_table. + call coupler_type_data_override('OCN', atm_fields, Time_import) + + ! Potentially override ice_ocean_boundary%fluxes from data_table. + ! Doing this before atmos_ocean_fluxes_calc call avoids unnecessary calculation of overridden fluxes. + ! However, we cannot use coupler_type_data_override here since it does not set the override flag on + ! overridden fields + do n = 1, ice_ocean_boundary%fluxes%num_bcs + do m = 1, ice_ocean_boundary%fluxes%bc(n)%num_fields + call data_override('OCN', ice_ocean_boundary%fluxes%bc(n)%field(m)%name, & + ice_ocean_boundary%fluxes%bc(n)%field(m)%values, Time_import, & + override=ice_ocean_boundary%fluxes%bc(n)%field(m)%override) + enddo + enddo + + ! Calculate the extra tracer fluxes + call get_domain_extent(ocean_public%domain, isc, iec, jsc, jec) + call atmos_ocean_fluxes_calc(atm_fields, ocean_public%fields, ice_ocean_boundary%fluxes, & + ice_ocean_boundary%ice_fraction, isc, iec, jsc, jec) + + ! Send diagnostics + call coupler_type_send_data(atm_fields, Time_import) + call coupler_type_send_data(ice_ocean_boundary%fluxes, Time_import) +#else call mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return +#endif !--------------- ! Update MOM6 @@ -1912,7 +2036,7 @@ subroutine ModelAdvance(gcomp, rc) ! determine restart filename call ESMF_ClockGetNextTime(clock, MyTime, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_TimeGet (MyTime, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc ) + call ESMF_TimeGet (MyTime, yy=year, mm=month, dd=day, s=seconds, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (cesm_coupled) then @@ -1975,6 +2099,14 @@ subroutine ModelAdvance(gcomp, rc) endif +#ifdef _USE_GENERIC_TRACER + ! Write fields for extra tracer fluxes to their internally defined ocean restart file + call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call gas_fields_restart(ocean_public%fields, ocean_public%domain, additional_restart_dir) +#endif + if (is_root_pe()) then write(stdout,*) subname//' writing restart file ',trim(restartname) endif @@ -2213,6 +2345,7 @@ subroutine ocean_model_finalize(gcomp, rc) type(ESMF_Alarm), allocatable :: alarmList(:) integer :: alarmCount logical :: write_restart + character(240) :: additional_restart_dir character(len=*),parameter :: subname='(MOM_cap:ocean_model_finalize)' real(8) :: MPI_Wtime, timefs @@ -2246,6 +2379,18 @@ subroutine ocean_model_finalize(gcomp, rc) call ocean_model_end(ocean_public, ocean_State, Time, write_restart=write_restart) +#ifdef _USE_GENERIC_TRACER + if (write_restart) then + ! Write fields for extra tracer fluxes to their internally defined ocean restart file + call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call gas_fields_restart(ocean_public%fields, ocean_public%domain, additional_restart_dir) + endif +#endif + + call field_manager_end() + call io_infra_end() call MOM_infra_end() diff --git a/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 b/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 new file mode 100644 index 0000000000..af1d886b37 --- /dev/null +++ b/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 @@ -0,0 +1,680 @@ +!> Contains routines for handling FMS coupler_bc_type tracer flux structures +!! when using MOM generic tracers + +module MOM_cap_gtracer_flux + +use MOM_domains, only: domain2d +use MOM_coupler_types, only: coupler_1d_bc_type, coupler_2d_bc_type +use MOM_coupler_types, only: ind_flux, ind_deltap, ind_kw, ind_flux0 +use MOM_coupler_types, only: ind_pcair, ind_u10, ind_psurf +use MOM_coupler_types, only: ind_alpha, ind_csurf, ind_sc_no +use MOM_coupler_types, only: ind_runoff, ind_deposition +use MOM_ocean_model_nuopc, only: ocean_model_flux_init +use coupler_types_mod, only: coupler_type_register_restarts, coupler_type_restore_state +use atmos_ocean_fluxes_mod, only: atmos_ocean_type_fluxes_init, atmos_ocean_fluxes_init +use fms2_io_mod, only: FmsNetcdfDomainFile_t +use fms2_io_mod, only: fms2_check_if_open => check_if_open, fms2_close_file => close_file +use fms2_io_mod, only: fms2_write_data => write_data +use fms2_io_mod, only: fms2_read_restart => read_restart, fms2_write_restart => write_restart +use fms2_io_mod, only: fms2_get_global_io_domain_indices => get_global_io_domain_indices +use field_manager_mod, only: fm_field_name_len, fm_type_name_len, fm_loop_over_list, fm_change_list +use fm_util_mod, only: fm_util_get_real_array +use mpp_mod, only: mpp_error, FATAL +use FMSconstants, only: wtmair, rdgas, vonkarm + +implicit none; private + +! Public member functions +public :: gas_exchange_init +public :: gas_fields_restore +public :: gas_fields_restart +public :: add_gas_fluxes_param +public :: get_coupled_field_name +public :: atmos_ocean_fluxes_calc +public :: UNKNOWN_CMEPS_FIELD + +character(len=*), parameter :: mod_name = 'mom_cap_gtracer_flux' +character(len=*), parameter :: UNKNOWN_CMEPS_FIELD = "UNKNOWN_FIELD" +real, parameter :: epsln=1.0e-30 + +!> FMS coupler_bc_types for additional tracer fields when using generic tracers +logical :: gas_fluxes_initialized = .false. ! This is set to true when the following types are initialized. +type(coupler_1d_bc_type), target :: ex_gas_fields_atm ! tracer fields in atm + !< Structure containing atmospheric surface variables that are used in the + !! calculation of the atmosphere-ocean gas fluxes, as well as parameters + !! regulating these fluxes. The fields in this structure are never actually + !! set, but the structure is used for initialisation of components and to + !! spawn other structure whose fields are set. +type(coupler_1d_bc_type), target :: ex_gas_fields_ocn ! tracer fields atop the ocean + !< Structure containing ocean surface variables that are used in the + !! calculation of the atmosphere-ocean gas fluxes, as well as parameters + !! regulating these fluxes. The fields in this structure are never actually + !! set, but the structure is used for initialisation of components and to + !! spawn other structure whose fields are set. +type(coupler_1d_bc_type), target :: ex_gas_fluxes ! tracer fluxes between the atm and ocean + !< A structure for exchanging gas or tracer fluxes between the atmosphere and + !! ocean, defined by the field table, as well as a place holder of + !! intermediate calculations, such as piston velocities, and parameters that + !! impact the fluxes. The fields in this structure are never actually set, + !! but the structure is used for initialisation of components and to spawn + !! other structure whose fields are set. + +contains + +!> \brief Gas and tracer field initialization routine for running with MOM generic tracers. +!! Copied and adapted slightly from +!! https://github.com/NOAA-GFDL/FMScoupler/blob/7761886/full/flux_exchange.F90#L626. +subroutine gas_exchange_init (gas_fields_atm, gas_fields_ocn, gas_fluxes) + type(coupler_1d_bc_type), optional, pointer :: gas_fields_atm ! tracer fields in atm + !< Pointer to a structure containing atmospheric surface variables that + !! are used in the calculation of the atmosphere-ocean gas fluxes, as well + !! as parameters regulating these fluxes. + type(coupler_1d_bc_type), optional, pointer :: gas_fields_ocn ! tracer fields atop the ocean + !< Pointer to a structure containing ocean surface variables that are + !! used in the calculation of the atmosphere-ocean gas fluxes, as well as + !! parameters regulating these fluxes. + type(coupler_1d_bc_type), optional, pointer :: gas_fluxes ! tracer fluxes between the atm and ocean + !< Pointer to a structure for exchanging gas or tracer fluxes between the + !! atmosphere and ocean, defined by the field table, as well as a place holder + !! of intermediate calculations, such as piston velocities, and parameters + !! that impact the fluxes. + + if (.not.gas_fluxes_initialized) then + call atmos_ocean_type_fluxes_init( ) + call ocean_model_flux_init( ) + call atmos_ocean_fluxes_init(ex_gas_fluxes, ex_gas_fields_atm, ex_gas_fields_ocn) + gas_fluxes_initialized = .true. + endif + + if (present(gas_fields_atm)) gas_fields_atm => ex_gas_fields_atm + if (present(gas_fields_ocn)) gas_fields_ocn => ex_gas_fields_ocn + if (present(gas_fluxes)) gas_fluxes => ex_gas_fluxes + +end subroutine gas_exchange_init + +!> \brief Restore FMS coupler_bc_type state from ocean restart file +! See https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1096 +subroutine gas_fields_restore(gas_fields, domain, directory) + type(coupler_2d_bc_type), intent(inout) :: gas_fields !< FMS coupler_bc_type to be registered for restarts + type(domain2D), intent(in) :: domain !< The FMS domain to use for this registration call + character(len=*), optional, intent(in) :: directory !< Directory containing the restart file + + ! local variables + type(FmsNetcdfDomainFile_t), dimension(:), pointer :: ocn_bc_restart => NULL() !< Structures describing + !! the restart files + integer :: num_ocn_bc_restart !< The number of restart + !! files to use + integer :: n + + call coupler_type_register_restarts(gas_fields, ocn_bc_restart, num_ocn_bc_restart, & + domain, to_read=.true., ocean_restart=.true., directory=directory) + + ! Restore the fields from the restart files + do n = 1, num_ocn_bc_restart + if (fms2_check_if_open(ocn_bc_restart(n))) then + call fms2_read_restart(ocn_bc_restart(n)) + endif + enddo + + ! Check whether the restarts were read successfully. + call coupler_type_restore_state(gas_fields, use_fms2_io=.true., test_by_field=.true.) + + do n = 1, num_ocn_bc_restart + if(fms2_check_if_open(ocn_bc_restart(n))) call fms2_close_file(ocn_bc_restart(n)) + enddo + +end subroutine gas_fields_restore + +!> \brief Write ocean restart file for FMS coupler_bc_type state +! See https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1408 +subroutine gas_fields_restart(gas_fields, domain, directory) + type(coupler_2d_bc_type), intent(inout) :: gas_fields !< FMS coupler_bc_type to be registered for restarts + type(domain2D), intent(in) :: domain !< The FMS domain to use for this registration call + character(len=*), optional, intent(in) :: directory !< Directory containing the restart file + + ! local variables + type(FmsNetcdfDomainFile_t), dimension(:), pointer :: ocn_bc_restart => NULL() !< Structures describing + !! the restart files + integer :: num_ocn_bc_restart !< The number of restart + !! files to use + integer :: n + + call coupler_type_register_restarts(gas_fields, ocn_bc_restart, num_ocn_bc_restart, & + domain, to_read=.false., ocean_restart=.true., directory=directory) + + do n = 1, num_ocn_bc_restart + if (fms2_check_if_open(ocn_bc_restart(n))) then + call fms2_write_restart(ocn_bc_restart(n)) + call add_domain_dimension_data(ocn_bc_restart(n)) + call fms2_close_file(ocn_bc_restart(n)) + endif + enddo + +end subroutine gas_fields_restart + +!> Register the axis data as a variable in the netcdf file and add some dummy data. +!! This is needed so the combiner can work correctly when the io_layout is not 1,1. Copied from +!! https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1328 +subroutine add_domain_dimension_data(fileobj) + type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2io domain decomposed fileobj + + ! local variables + integer, dimension(:), allocatable :: buffer !< Buffer with axis data + integer :: is, ie !< Starting and Ending indices for data + + call fms2_get_global_io_domain_indices(fileobj, "xaxis_1", is, ie, indices=buffer) + call fms2_write_data(fileobj, "xaxis_1", buffer) + deallocate(buffer) + + call fms2_get_global_io_domain_indices(fileobj, "yaxis_1", is, ie, indices=buffer) + call fms2_write_data(fileobj, "yaxis_1", buffer) + deallocate(buffer) + +end subroutine add_domain_dimension_data + +!> Retrieve param array from field_manager and add to FMS coupler_bc_type. This is +!! needed because the coupler_type_spawn routine does not copy the param array into +!! the spawned type. Hopefully we can get rid of this. This routine is based on +!! https://github.com/NOAA-GFDL/FMS/blob/7f58528/coupler/atmos_ocean_fluxes.F90#L448 +subroutine add_gas_fluxes_param(gas_fluxes) + type(coupler_2d_bc_type), intent(inout) :: gas_fluxes !< FMS coupler_bc_type to add param to + + ! local variables + integer :: n + character(len=fm_field_name_len) :: name + character(len=fm_type_name_len) :: typ + integer :: ind + + character(len=*), parameter :: sub_name = 'add_gas_fluxes_param' + character(len=*), parameter :: error_header =& + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + + n = 0 + do while (fm_loop_over_list('/coupler_mod/fluxes', name, typ, ind)) + if (typ .ne. 'list') then + call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list') + endif + + n = n + 1 + + if (.not. fm_change_list('/coupler_mod/fluxes/' // trim(name))) then + call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name)) + endif + + if (gas_fluxes%bc(n)%name .eq. name) then + gas_fluxes%bc(n)%param => fm_util_get_real_array('param') + else + call mpp_error(FATAL, trim(error_header) // ' Problem setting param array pointer') + endif + enddo +end subroutine add_gas_fluxes_param + +!> Return the CMEPS standard_name of the coupled field required for a given coupled +!! generic_tracer flux name. +function get_coupled_field_name(name) + character(len=64) :: get_coupled_field_name !< CMEPS standard_name + character(len=*), intent(in) :: name !< gtracer flux name + + ! Add other coupled field names here + select case(trim(name)) + case( 'co2_flux' ) + get_coupled_field_name = "Sa_co2prog" + case default + get_coupled_field_name = UNKNOWN_CMEPS_FIELD + end select +end function get_coupled_field_name + +!> \brief Calculate the FMS coupler_bc_type ocean tracer fluxes. Units should be mol/m^2/s. +!! Upward flux is positive. +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +!! and subsequently modified in the following ways: +!! - Operate on 2D inputs, rather than 1D +!! - Add calculation for 'air_sea_deposition' taken from +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_dep_fluxes_calc.F90 +!! - Multiply fluxes by ice_fraction input, rather than masking based on seawater input +!! - Use MOM over FMS modules where easy to do so +!! - Make tsurf input optional, as it is only used by a few implementations +!! - Use ind_runoff rather than ind_deposition in runoff flux calculation (note, their +!! values are equal) +!! - Rename gas_fields_ice to gas_fields_ocn +subroutine atmos_ocean_fluxes_calc(gas_fields_atm, gas_fields_ocn, gas_fluxes,& + ice_fraction, isc, iec, jsc, jec, tsurf, ustar, cd_m) + type(coupler_2d_bc_type), intent(in) :: gas_fields_atm ! fields in atm + !< Structure containing atmospheric surface variables that are used in the calculation + !! of the atmosphere-ocean tracer fluxes. + type(coupler_2d_bc_type), intent(in) :: gas_fields_ocn ! fields atop the ocean + !< Structure containing ocean surface variables that are used in the calculation of the + !! atmosphere-ocean tracer fluxes. + type(coupler_2d_bc_type), intent(inout) :: gas_fluxes ! fluxes between the atm and ocean + !< Structure containing the gas fluxes between the atmosphere and the ocean and + !! parameters related to the calculation of these fluxes. + real, intent(in) :: ice_fraction(isc:iec,jsc:jec) !< sea ice fraction + integer, intent(in) :: isc !< The start i-index of cell centers within + !! the computational domain + integer, intent(in) :: iec !< The end i-index of cell centers within the + !! computational domain + integer, intent(in) :: jsc !< The start j-index of cell centers within + !! the computational domain + integer, intent(in) :: jec !< The end j-index of cell centers within the + !! computational domain + real, intent(in), optional :: tsurf(isc:iec,jsc:jec) !< surface temperature + real, intent(in), optional :: ustar(isc:iec,jsc:jec) !< friction velocity, not + !! used + real, intent(in), optional :: cd_m (isc:iec,jsc:jec) !< drag coefficient, not + !! used + + ! local variables + character(len=*), parameter :: sub_name = 'atmos_ocean_fluxes_calc' + character(len=*), parameter :: error_header =& + & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + real, parameter :: permeg=1.0e-6 + + integer :: n + integer :: i + integer :: j + real, dimension(:,:), allocatable :: kw + real, dimension(:,:), allocatable :: cair + character(len=128) :: error_string + + ! Return if no fluxes to be calculated + if (gas_fluxes%num_bcs .le. 0) return + + if (.not. associated(gas_fluxes%bc)) then + if (gas_fluxes%num_bcs .ne. 0) then + call mpp_error(FATAL, trim(error_header) // ' Number of gas fluxes not zero') + else + return + endif + endif + + do n = 1, gas_fluxes%num_bcs + ! only do calculations if the flux has not been overridden + if ( .not. gas_fluxes%bc(n)%field(ind_flux)%override) then + if (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux_generic') then + if (.not. allocated(kw)) then + allocate( kw(isc:iec,jsc:jec) ) + allocate ( cair(isc:iec,jsc:jec) ) + elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then + call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match') + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =& + & (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) * & + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2 + cair(i,j) = & + gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & sqrt(660. / (gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j) + epsln)) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & sqrt(660. / (gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j) + epsln)) *& + & gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) / & + (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'duce') then + if (.not. present(tsurf)) then + call mpp_error(FATAL, trim(error_header) // ' Implementation ' //& + trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //& + ' requires input tsurf') + endif + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) = & + & (1 - ice_fraction(i,j)) * gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) /& + & (770.+45.*gas_fluxes%bc(n)%param(1)**(1./3.)) *& + & 101325./(rdgas*wtmair*1e-3*tsurf(i,j) *& + & max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln)) + !alpha: mol/m3/atm + cair(i,j) = & + gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6 + cair(i,j) = max(cair(i,j),0.) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /& + & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'johnson') then + if (.not. present(tsurf)) then + call mpp_error(FATAL, trim(error_header) // ' Implementation ' //& + trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //& + ' requires input tsurf') + endif + !f1p: not sure how to pass salinity. For now, just force at 35. + do j = jsc,jec + do i = isc,iec + !calc_kw(tk,p,u10,h,vb,mw,sc_w,ustar,cd_m) + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =& + & (1 - ice_fraction(i,j)) * calc_kw(tsurf(i,j),& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j),& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j),& + & 101325./(rdgas*wtmair*1e-3*tsurf(i,j)*& + & max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln)),& + & gas_fluxes%bc(n)%param(2),& + & gas_fluxes%bc(n)%param(1),& + & gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j)) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6 + cair(i,j) = max(cair(i,j),0.) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /& + & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux') then + if (.not. allocated(kw)) then + allocate( kw(isc:iec,jsc:jec) ) + allocate ( cair(isc:iec,jsc:jec) ) + elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then + call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match') + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2_data') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2 + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'linear') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & max(0.0, gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) - gas_fluxes%bc(n)%param(2)) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(3) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_deposition') then + if (gas_fluxes%bc(n)%param(1) .le. 0.0) then + write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1) + call mpp_error(FATAL, 'Bad parameter (' // trim(error_string) //& + & ') for air_sea_deposition for ' // trim(gas_fluxes%bc(n)%name)) + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'dry') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'wet') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1) + enddo + enddo + else + call mpp_error(FATAL, 'Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'land_sea_runoff') then + if (gas_fluxes%bc(n)%param(1) .le. 0.0) then + write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1) + call mpp_error(FATAL, ' Bad parameter (' // trim(error_string) //& + & ') for land_sea_runoff for ' // trim(gas_fluxes%bc(n)%name)) + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'river') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + & gas_fields_atm%bc(n)%field(ind_runoff)%values(i,j) /& + & gas_fluxes%bc(n)%param(1) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + else + call mpp_error(FATAL, ' Unknown flux_type (' // trim(gas_fluxes%bc(n)%flux_type) //& + & ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + endif + enddo + + if (allocated(kw)) then + deallocate(kw) + deallocate(cair) + endif +end subroutine atmos_ocean_fluxes_calc + +!> Calculate \f$k_w\f$ +!! +!! Taken from Johnson, Ocean Science, 2010. (http://doi.org/10.5194/os-6-913-2010) +!! +!! Uses equations defined in Liss[1974], +!! \f[ +!! F = K_g(c_g - H C_l) = K_l(c_g/H - C_l) +!! \f] +!! where \f$c_g\f$ and \f$C_l\f$ are the bulk gas and liquid concentrations, \f$H\f$ +!! is the Henry's law constant (\f$H = c_{sg}/C_{sl}\f$, where \f$c_{sg}\f$ is the +!! equilibrium concentration in gas phase (\f$g/cm^3\f$ of air) and \f$C_{sl}\f$ is the +!! equilibrium concentration of unionised dissolved gas in liquid phase (\f$g/cm^3\f$ +!! of water)), +!! \f[ +!! 1/K_g = 1/k_g + H/k_l +!! \f] +!! and +!! \f[ +!! 1/K_l = 1/k_l + 1/{Hk_g} +!! \f] +!! where \f$k_g\f$ and \f$k_l\f$ are the exchange constants for the gas and liquid +!! phases, respectively. +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function calc_kw(tk, p, u10, h, vb, mw, sc_w, ustar, cd_m) + real, intent(in) :: tk !< temperature at surface in kelvin + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s + real, intent(in) :: h !< Henry's law constant (\f$H=c_sg/C_sl\f$) (unitless) + real, intent(in) :: vb !< Molar volume + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: sc_w + real, intent(in), optional :: ustar !< Friction velocity (m/s). If not provided, + !! ustar = \f$u_{10} \sqrt{C_D}\f$. + real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$). Used only if + !! ustar is provided. + !! If ustar is not provided, + !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$ + + real :: ra,rl,tc + + tc = tk-273.15 + ra = 1./max(h*calc_ka(tc,p,mw,vb,u10,ustar,cd_m),epsln) + rl = 1./max(calc_kl(tc,u10,sc_w),epsln) + calc_kw = 1./max(ra+rl,epsln) +end function calc_kw + +!> Calculate \f$k_a\f$ +!! +!! See calc_kw +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function calc_ka(t, p, mw, vb, u10, ustar, cd_m) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< molar volume + real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s + real, intent(in), optional :: ustar !< Friction velocity (m/s). If not provided, + !! ustar = \f$u_{10} \sqrt{C_D}\f$. + real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$). Used only if + !! ustar is provided. + !! If ustar is not provided, + !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$ + + real :: sc + real :: ustar_t, cd_m_t + + if (.not. present(ustar)) then + !drag coefficient + cd_m_t = 6.1e-4 +0.63e-4*u10 + !friction velocity + ustar_t = u10*sqrt(cd_m_t) + else + cd_m_t = cd_m + ustar_t = ustar + end if + sc = schmidt_g(t,p,mw,vb) + calc_ka = 1e-3+ustar_t/(13.3*sqrt(sc)+1/sqrt(cd_m_t)-5.+log(sc)/(2.*vonkarm)) +end function calc_ka + +!> Calculate \f$k_l\f$ +!! +!! See calc_kw, and Nightingale, Global Biogeochemical Cycles, 2000 +!! (https://doi.org/10.1029/1999GB900091) +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function calc_kl(t, v, sc) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: v !< wind speed at surface in m/s + real, intent(in) :: sc + + calc_kl = (((0.222*v**2)+0.333*v)*(max(sc,epsln)/600.)**(-0.5))/(100.*3600.) +end function calc_kl + +!> Schmidt number of the gas in air +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function schmidt_g(t, p, mw, vb) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< molar volume + + real :: d,v + + d = d_air(t,p,mw,vb) + v = v_air(t) + schmidt_g = v / d +end function schmidt_g + +!> From Fuller, Industrial & Engineering Chemistry (https://doi.org/10.1021/ie50677a007) +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function d_air(t, p, mw, vb) + real, intent(in) :: t !< temperature in c + real, intent(in) :: p !< pressure in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< diffusion coefficient (\f$cm3/mol\f$) + + real, parameter :: ma = 28.97d0 !< molecular weight air in g/mol + real, parameter :: va = 20.1d0 !< diffusion volume for air (\f$cm^3/mol\f$) + + real :: pa + + ! convert p to atm + pa = 9.8692d-6*p + d_air = 1d-3 *& + & (t+273.15d0)**(1.75d0)*sqrt(1d0/ma + 1d0/mw)/(pa*(va**(1d0/3d0)+vb**(1d0/3d0))**2d0) + ! d_air is in cm2/s convert to m2/s + d_air = d_air * 1d-4 +end function d_air + +!> kinematic viscosity in air +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function p_air(t) + real, intent(in) :: t + + real, parameter :: sd_0 = 1.293393662d0,& + & sd_1 = -5.538444326d-3,& + & sd_2 = 3.860201577d-5,& + & sd_3 = -5.2536065d-7 + p_air = sd_0+(sd_1*t)+(sd_2*t**2)+(sd_3*t**3) +end function p_air + +!> Kinematic viscosity in air (\f$m^2/s\f$ +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function v_air(t) + real, intent(in) :: t !< temperature in C + v_air = n_air(t)/p_air(t) +end function v_air + +!> dynamic viscosity in air +!! +!! This routine was copied from FMScoupler at +!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90 +real function n_air(t) + real, intent(in) :: t !< temperature in C + + real, parameter :: sv_0 = 1.715747771d-5,& + & sv_1 = 4.722402075d-8,& + & sv_2 = -3.663027156d-10,& + & sv_3 = 1.873236686d-12,& + & sv_4 = -8.050218737d-14 + ! in n.s/m^2 (pa.s) + n_air = sv_0+(sv_1*t)+(sv_2*t**2)+(sv_3*t**3)+(sv_4*t**4) +end function n_air + +end module MOM_cap_gtracer_flux \ No newline at end of file diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 809a507e5e..cde8727499 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -20,8 +20,15 @@ module MOM_cap_methods use MOM_surface_forcing_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type use MOM_domains, only: pass_var +use MOM_coupler_types, only: coupler_2d_bc_type use mpp_domains_mod, only: mpp_get_compute_domain +#ifdef _USE_GENERIC_TRACER +use MOM_coupler_types, only: set_coupler_type_data +use MOM_coupler_types, only: ind_pcair, ind_u10, ind_psurf, ind_runoff, ind_deposition +use MOM_cap_gtracer_flux, only: get_coupled_field_name, UNKNOWN_CMEPS_FIELD +#endif + ! By default make data private implicit none; private @@ -72,11 +79,17 @@ end subroutine mom_set_geomtype !> This function has a few purposes: !! (1) it imports surface fluxes using data from the mediator; and !! (2) it can apply restoring in SST and SSS. -subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, rc) +!! (3) optional: if atm_fields is provided, it imports and sets the fields in atm_fields required +!! for the calculation of coupled generic tracer fluxes +subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, atm_fields, rc) type(ocean_public_type) , intent(in) :: ocean_public !< Ocean surface state type(ocean_grid_type) , intent(in) :: ocean_grid !< Ocean model grid type(ESMF_State) , intent(inout) :: importState !< incoming data from mediator type(ice_ocean_boundary_type) , intent(inout) :: ice_ocean_boundary !< Ocean boundary forcing + type(coupler_2d_bc_type), optional, intent(inout) :: atm_fields !< If present, this type + !! describes the atmospheric tracer fields to + !! be imported for the calculation of generic + !! tracer fluxes. integer , intent(inout) :: rc !< Return code ! Local Variables @@ -89,7 +102,10 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, real(ESMF_KIND_R8), allocatable :: tauy(:,:) real(ESMF_KIND_R8), allocatable :: stkx(:,:,:) real(ESMF_KIND_R8), allocatable :: stky(:,:,:) + real(ESMF_KIND_R8), allocatable :: work(:,:) character(len=*) , parameter :: subname = '(mom_import)' + character(len=256) :: stdname + integer :: field_index rc = ESMF_SUCCESS @@ -551,6 +567,48 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, deallocate(stkx,stky) endif + !--- + ! Tracer flux fields for generic tracers + !--- +#ifdef _USE_GENERIC_TRACER + if (present(atm_fields)) then + ! Set fields in atm_fields from coupler + allocate (work(isc:iec,jsc:jec)) + do n = 1, atm_fields%num_bcs + if (atm_fields%bc(n)%flux_type .eq. 'air_sea_deposition') then + field_index = ind_deposition + elseif (atm_fields%bc(n)%flux_type .eq. 'land_sea_runoff') then + field_index = ind_runoff + else + ! This is a gas flux - set ind_u10 and ind_psurf + ! Note, we set these fields even though the pcair field may not be set below. This + ! is to allow flux calculation with overridden pcair fields + field_index = ind_pcair + call set_coupler_type_data(sqrt(ice_ocean_boundary%u10_sqr), n, atm_fields, & + idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=ind_u10) + call set_coupler_type_data(ice_ocean_boundary%p, n, atm_fields, & + idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=ind_psurf) + endif + + stdname = get_coupled_field_name(atm_fields%bc(n)%name) + if (stdname /= UNKNOWN_CMEPS_FIELD) then + call ESMF_LogWrite(trim(subname)//': generic_tracer flux, '//trim(atm_fields%bc(n)%name)//& + ': setting field index '//CHAR(48+field_index)//' to '//trim(stdname)//' if provided '//& + 'by coupler, otherwise defaulting to zero', ESMF_LOGMSG_INFO) + work(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, trim(stdname), isc, iec, jsc, jec, work, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call set_coupler_type_data(work, n, atm_fields, & + idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=field_index) + else + call ESMF_LogWrite(trim(subname)//': generic_tracer flux, '//trim(atm_fields%bc(n)%name)//& + ': no fields set from coupler', ESMF_LOGMSG_INFO) + endif + enddo + deallocate(work) + endif +#endif + end subroutine mom_import !> Maps outgoing ocean data to ESMF State @@ -615,7 +673,8 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, jg = j + ocean_grid%jsc - jsc do i = isc, iec ig = i + ocean_grid%isc - isc - omask(i,j) = nint(ocean_grid%mask2dT(ig,jg)) + !omask(i,j) = nint(ocean_grid%mask2dT(ig,jg)) - ceiling(ocean_state%ice_shelf_CSp%ISS%hmask(ig,jg)) + omask(i,j) = nint(ocean_grid%mask2dT(ig,jg)) - ceiling(ocean_state%fluxes%frac_shelf_h(ig,jg)) enddo enddo diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index a83576028a..b46cbaa029 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -20,6 +20,7 @@ module MOM_ocean_model_nuopc use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf use MOM_diag_mediator, only : diag_ctrl, enable_averages, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end +use MOM_domains, only : MOM_domain_type, domain2d, clone_MOM_domain, get_domain_extent use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE use MOM_domains, only : TO_ALL, Omit_Corners use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe @@ -46,13 +47,14 @@ module MOM_ocean_model_nuopc use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use MOM_ice_shelf, only : ice_shelf_query +use MOM_data_override, only : data_override_init use MOM_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type use MOM_coupler_types, only : coupler_type_spawn, coupler_type_write_chksums use MOM_coupler_types, only : coupler_type_initialized, coupler_type_copy_data use MOM_coupler_types, only : coupler_type_set_diags, coupler_type_send_data -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain -use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use MOM_io, only : stdout use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init @@ -102,7 +104,7 @@ module MOM_ocean_model_nuopc !! points of the two velocity components. Valid entries !! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, !! corresponding to the community-standard Arakawa notation. - !! (These are named integers taken from mpp_parameter_mod.) + !! (These are named integers taken from the MOM_domains module.) !! Following MOM5, stagger is BGRID_NE by default when the !! ocean is initialized, but here it is set to -999 so that !! a global max across ocean and non-ocean processors can be @@ -133,7 +135,7 @@ module MOM_ocean_model_nuopc !> The ocean_state_type contains all information about the state of the ocean, !! with a format that is private so it can be readily changed without disrupting !! other coupled components. -type, public :: ocean_state_type ; private +type, public :: ocean_state_type ! This type is private, and can therefore vary between different ocean models. logical :: is_ocean_PE = .false. !< True if this is an ocean PE. type(time_type) :: Time !< The ocean model's time and master clock. @@ -209,6 +211,7 @@ module MOM_ocean_model_nuopc Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This !! is null if there is no ice shelf. + logical :: override_melt !< If true, override melt using data_override type(marine_ice_CS), pointer :: & marine_ice_CSp => NULL() !< A pointer to the control structure for the !! marine ice effects module. @@ -387,15 +390,24 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! vertical integrals, since the related 3-d sums are not negligible in cost. call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, & do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & - use_meltpot=use_melt_pot, use_marbl_tracers=OS%use_MARBL) + use_meltpot=use_melt_pot, use_iceshelves=OS%use_ice_shelf, & + use_marbl_tracers=OS%use_MARBL) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves) if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & - OS%diag, Time_init, OS%dirs%output_directory, OS%forces, OS%fluxes) + OS%diag, Time_init, OS%dirs%output_directory) !, OS%forces, OS%fluxes) endif + if (OS%use_ice_shelf) then + call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes) + call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%flux_tmp) + call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces) + call ice_shelf_query(OS%ice_shelf_CSp, OS%grid, data_override_melt=OS%override_melt) + if (OS%override_melt) call data_override_init(Ocean_Domain_in=OS%grid%domain%mpp_domain) + endif + if (OS%icebergs_alter_ocean) then call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp) if (.not. OS%use_ice_shelf) & @@ -411,14 +423,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! it also initializes statistical waves. call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) - if (associated(OS%grid%Domain%maskmap)) then - call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & - OS%diag, maskmap=OS%grid%Domain%maskmap, & - gas_fields_ocn=gas_fields_ocn) - else - call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & - OS%diag, gas_fields_ocn=gas_fields_ocn) - endif + call initialize_ocean_public_type(OS%grid%Domain, Ocean_sfc, OS%diag, gas_fields_ocn=gas_fields_ocn) ! This call can only occur here if the coupler_bc_type variables have been ! initialized already using the information from gas_fields_ocn. @@ -533,8 +538,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) ! Translate Ice_ocean_boundary into fluxes. - call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), & - index_bnds(3), index_bnds(4)) + call get_domain_extent(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), index_bnds(3), index_bnds(4)) weight = 1.0 @@ -747,7 +751,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu OS%dirs%restart_output_dir) ! Is this needed? if (OS%use_ice_shelf) then call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, & - OS%dirs%restart_output_dir) + './RESTART/') !OS%dirs%restart_output_dir) endif else if (BTEST(OS%Restart_control,1)) then @@ -756,7 +760,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, time_stamped=.true.) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir, .true.) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/', .true.) ! OS%dirs%restart_output_dir, .true.) endif endif if (BTEST(OS%Restart_control,0)) then @@ -765,7 +769,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/') !OS%dirs%restart_output_dir) endif endif endif @@ -800,7 +804,7 @@ end subroutine ocean_model_end subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) type(ocean_state_type), pointer :: OS !< A pointer to the structure containing the !! internal ocean state (in). - type(time_type), intent(in) :: Time !< The model time at this call, needed for mpp_write calls. + type(time_type), intent(in) :: Time !< The model time at this call, needed for writing files. character(len=*), optional, intent(in) :: directory !< An optional directory into which to !! write these restart files. character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-stamp) @@ -826,21 +830,17 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/') ! OS%dirs%restart_output_dir) endif end subroutine ocean_model_save_restart !> Initialize the public ocean type -subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, & - gas_fields_ocn) - type(domain2D), intent(in) :: input_domain !< The ocean model domain description +subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_fields_ocn) + type(MOM_domain_type), intent(in) :: input_domain !< The ocean model domain description type(ocean_public_type), intent(inout) :: Ocean_sfc !< A structure containing various publicly - !! visible ocean surface properties after initialization, whose - !! elements are allocated here. - type(diag_ctrl), intent(in) :: diag !< A structure that regulates diagnsotic output - logical, dimension(:,:), & - optional, intent(in) :: maskmap !< A mask indicating which virtual processors - !! are actually in use. If missing, all are used. + !! visible ocean surface properties after + !! initialization, whose elements are allocated here. + type(diag_ctrl), intent(in) :: diag !< A structure that regulates diagnostic output type(coupler_1d_bc_type), & optional, intent(in) :: gas_fields_ocn !< If present, this type describes the !! ocean and surface-ice fields that will participate @@ -852,14 +852,9 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, ! and have no halos. integer :: isc, iec, jsc, jec - call mpp_get_layout(input_domain,layout) - call mpp_get_global_domain(input_domain, xsize=xsz, ysize=ysz) - if (PRESENT(maskmap)) then - call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain, maskmap=maskmap) - else - call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain) - endif - call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec) + call clone_MOM_domain(input_domain, Ocean_sfc%Domain, halo_size=0, symmetric=.false.) + + call get_domain_extent(Ocean_sfc%Domain, isc, iec, jsc, jec) allocate(Ocean_sfc%t_surf (isc:iec,jsc:jec), & ! time averaged sst (Kelvin) passed to atmosphere/ice model Ocean_sfc%s_surf (isc:iec,jsc:jec), & ! time averaged sss (psu) passed to atmosphere/ice models @@ -909,8 +904,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec call pass_vector(sfc_state%u, sfc_state%v, G%Domain) - call mpp_get_compute_domain(Ocean_sfc%Domain, isc_bnd, iec_bnd, & - jsc_bnd, jec_bnd) + call get_domain_extent(Ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd) if (present(patm)) then ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd). if (.not.present(press_to_z)) call MOM_error(FATAL, & diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 index 42c386497a..90e7d795ff 100644 --- a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 +++ b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 @@ -21,6 +21,7 @@ module generic_tracer public generic_tracer_vertdiff_G public generic_tracer_get_diag_list public generic_tracer_coupler_accumulate + public generic_tracer_update_from_coupler !> Turn on generic tracers (note dangerous use of module data) logical :: do_generic_tracer = .true. @@ -65,6 +66,14 @@ subroutine generic_tracer_coupler_accumulate(IOB_struc, weight, model_time) type(time_type), optional,intent(in) :: model_time !< Time end subroutine generic_tracer_coupler_accumulate + !> Modify the values obtained from the coupler + subroutine generic_tracer_update_from_coupler(ilb, jlb, salt_flux_added) + integer, intent(in) :: ilb !< Lower bounds of x extent of input arrays on data domain + integer, intent(in) :: jlb !< Lower bounds of y extent of input arrays on data domain + real, dimension(ilb:,jlb:), intent(in) :: salt_flux_added !< Surface salt flux into ocean from restoring + !! or flux adjustment [g/m^2/sec] + end subroutine generic_tracer_update_from_coupler + !> Calls the corresponding generic_X_update_from_source routine for each package X subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dtts,& grid_dat,model_time,nbands,max_wavelength_band,sw_pen_band,opacity_band,internal_heat,& diff --git a/config_src/infra/FMS1/MOM_couplertype_infra.F90 b/config_src/infra/FMS1/MOM_couplertype_infra.F90 index 637f2b5ebf..0fb57c991a 100644 --- a/config_src/infra/FMS1/MOM_couplertype_infra.F90 +++ b/config_src/infra/FMS1/MOM_couplertype_infra.F90 @@ -9,7 +9,10 @@ module MOM_couplertype_infra use coupler_types_mod, only : coupler_type_increment_data, coupler_type_rescale_data use coupler_types_mod, only : coupler_type_extract_data, coupler_type_set_data use coupler_types_mod, only : coupler_type_data_override -use coupler_types_mod, only : ind_flux, ind_alpha, ind_csurf +use coupler_types_mod, only : ind_flux, ind_deltap, ind_kw, ind_flux0 +use coupler_types_mod, only : ind_pcair, ind_u10, ind_psurf +use coupler_types_mod, only : ind_alpha, ind_csurf, ind_sc_no +use coupler_types_mod, only : ind_runoff, ind_deposition use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux use MOM_domain_infra, only : domain2D @@ -22,7 +25,10 @@ module MOM_couplertype_infra public :: CT_set_data, CT_increment_data, CT_rescale_data public :: CT_copy_data, CT_extract_data, CT_redistribute_data public :: atmos_ocn_coupler_flux -public :: ind_flux, ind_alpha, ind_csurf +public :: ind_flux, ind_deltap, ind_kw, ind_flux0 +public :: ind_pcair, ind_u10, ind_psurf +public :: ind_alpha, ind_csurf, ind_sc_no +public :: ind_runoff, ind_deposition public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type !> This is the interface to spawn one coupler_bc_type into another. diff --git a/config_src/infra/FMS2/MOM_couplertype_infra.F90 b/config_src/infra/FMS2/MOM_couplertype_infra.F90 index 3bcccc1dc7..50fd5b09cd 100644 --- a/config_src/infra/FMS2/MOM_couplertype_infra.F90 +++ b/config_src/infra/FMS2/MOM_couplertype_infra.F90 @@ -9,7 +9,10 @@ module MOM_couplertype_infra use coupler_types_mod, only : coupler_type_increment_data, coupler_type_rescale_data use coupler_types_mod, only : coupler_type_extract_data, coupler_type_set_data use coupler_types_mod, only : coupler_type_data_override -use coupler_types_mod, only : ind_flux, ind_alpha, ind_csurf +use coupler_types_mod, only : ind_flux, ind_deltap, ind_kw, ind_flux0 +use coupler_types_mod, only : ind_pcair, ind_u10, ind_psurf +use coupler_types_mod, only : ind_alpha, ind_csurf, ind_sc_no +use coupler_types_mod, only : ind_runoff, ind_deposition use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux use MOM_domain_infra, only : domain2D @@ -22,7 +25,10 @@ module MOM_couplertype_infra public :: CT_set_data, CT_increment_data, CT_rescale_data public :: CT_copy_data, CT_extract_data, CT_redistribute_data public :: atmos_ocn_coupler_flux -public :: ind_flux, ind_alpha, ind_csurf +public :: ind_flux, ind_deltap, ind_kw, ind_flux0 +public :: ind_pcair, ind_u10, ind_psurf +public :: ind_alpha, ind_csurf, ind_sc_no +public :: ind_runoff, ind_deposition public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type !> This is the interface to spawn one coupler_bc_type into another. diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 156a397ff6..77d24a5ecb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -4077,7 +4077,7 @@ subroutine extract_surface_state(CS, sfc_state_in) depth_ml = CS%Hmix_UV if (CS%answer_date < 20190101) depth_ml = GV%H_to_Z*CS%Hmix_UV !$OMP parallel do default(shared) private(depth,dh,hv) - do J=js-1,ie + do J=js-1,je do i=is,ie depth(i) = 0.0 sfc_state%v(i,J) = 0.0 diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index f91d958fe8..2c3744cfa0 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -2421,6 +2421,11 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces) fluxes%salt_flux(i,j) = wt1*fluxes%salt_flux(i,j) + wt2*flux_tmp%salt_flux(i,j) enddo ; enddo + if (associated(fluxes%salt_flux_added) .and. associated(flux_tmp%salt_flux_added)) then + do j=js,je ; do i=is,ie + fluxes%salt_flux_added(i,j) = wt1*fluxes%salt_flux_added(i,j) + wt2*flux_tmp%salt_flux_added(i,j) + enddo ; enddo + endif if (associated(fluxes%heat_added) .and. associated(flux_tmp%heat_added)) then do j=js,je ; do i=is,ie fluxes%heat_added(i,j) = wt1*fluxes%heat_added(i,j) + wt2*flux_tmp%heat_added(i,j) diff --git a/src/framework/MOM_coupler_types.F90 b/src/framework/MOM_coupler_types.F90 index 25a2937aaa..82d273949a 100644 --- a/src/framework/MOM_coupler_types.F90 +++ b/src/framework/MOM_coupler_types.F90 @@ -9,7 +9,10 @@ module MOM_coupler_types use MOM_couplertype_infra, only : CT_copy_data, CT_increment_data, CT_rescale_data use MOM_couplertype_infra, only : CT_set_data, CT_extract_data, CT_redistribute_data use MOM_couplertype_infra, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type -use MOM_couplertype_infra, only : ind_flux, ind_alpha, ind_csurf +use MOM_couplertype_infra, only : ind_flux, ind_deltap, ind_kw, ind_flux0 +use MOM_couplertype_infra, only : ind_pcair, ind_u10, ind_psurf +use MOM_couplertype_infra, only : ind_alpha, ind_csurf, ind_sc_no +use MOM_couplertype_infra, only : ind_runoff, ind_deposition use MOM_domain_infra, only : domain2D use MOM_time_manager, only : time_type @@ -23,7 +26,10 @@ module MOM_coupler_types public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type ! These are encoding constant parameters that indicate whether a flux, solubility or ! surface ocean concentration are being set or accessed with an inquiry. -public :: ind_flux, ind_alpha, ind_csurf +public :: ind_flux, ind_deltap, ind_kw, ind_flux0 +public :: ind_pcair, ind_u10, ind_psurf +public :: ind_alpha, ind_csurf, ind_sc_no +public :: ind_runoff, ind_deposition !> This is the interface to spawn one coupler_bc_type into another. interface coupler_type_spawn diff --git a/src/framework/MOM_netcdf.F90 b/src/framework/MOM_netcdf.F90 index 4a7a61ec1c..ccdc399f5f 100644 --- a/src/framework/MOM_netcdf.F90 +++ b/src/framework/MOM_netcdf.F90 @@ -19,6 +19,7 @@ module MOM_netcdf use netcdf, only : nf90_inquire, nf90_inquire_dimension, nf90_inquire_variable use netcdf, only : nf90_inq_dimids, nf90_inq_varids use netcdf, only : NF90_MAX_NAME +use netcdf, only : nf90_inq_varid use MOM_error_handler, only : MOM_error, FATAL use MOM_io_infra, only : READONLY_FILE, WRITEONLY_FILE diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index b99cd3f184..acc27f0bd6 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -2109,7 +2109,7 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, integer :: nf ! The number of files that have been found so far integer :: m, length logical :: still_looking ! If true, the code is still looking for automatically named files - logical :: fexists ! True if a file has been found + logical :: fexists, fexists_decomp ! True if a file has been found character(len=32) :: filename_appendix = '' ! Filename appendix for ensemble runs character(len=80) :: restartname @@ -2190,9 +2190,17 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, else filepath = trim(directory)//trim(fname) inquire(file=filepath, exist=fexists) - if (.not. fexists) filepath = trim(filepath)//".nc" + if (.not.fexists .and. CS%parallel_restartfiles) & + fexists_decomp = file_exists(filepath, G%Domain) + + ! If not found, try with ".nc" extension + if (.not.(fexists .or. fexists_decomp)) then + filepath = trim(filepath)//".nc" + inquire(file=filepath, exist=fexists) + if (.not.fexists .and. CS%parallel_restartfiles) & + fexists_decomp = file_exists(filepath, G%Domain) + endif - inquire(file=filepath, exist=fexists) if (fexists) then nf = nf + 1 if (present(IO_handles)) & @@ -2202,6 +2210,14 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, if (present(file_paths)) file_paths(nf) = filepath if (is_root_pe() .and. (present(IO_handles))) & call MOM_error(NOTE,"MOM_restart: MOM run restarted using : "//trim(filepath)) + elseif (fexists_decomp) then + nf = nf + 1 + if (present(IO_handles)) & + call IO_handles(nf)%open(trim(filepath), READONLY_FILE, MOM_domain=G%Domain) + if (present(global_files)) global_files(nf) = .false. + if (present(file_paths)) file_paths(nf) = filepath + if (is_root_pe() .and. present(IO_handles)) & + call MOM_error(NOTE, "MOM_restart: MOM run restarted using decomposed fileset: "//trim(filepath)) else if (present(IO_handles)) & call MOM_error(WARNING,"MOM_restart: Unable to find restart file : "//trim(filepath)) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 2def8097ea..90d81284a1 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -86,7 +86,7 @@ module MOM_ice_shelf ! vary with the Boussinesq approximation, the Boussinesq variant is given first. !> Control structure that contains ice shelf parameters and diagnostics handles -type, public :: ice_shelf_CS ; private +type, public :: ice_shelf_CS ! Parameters type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control !! structure for the ice shelves @@ -126,9 +126,9 @@ module MOM_ice_shelf real :: kd_molec_salt!< The molecular diffusivity of salt [Z2 T-1 ~> m2 s-1]. real :: kd_molec_temp!< The molecular diffusivity of heat [Z2 T-1 ~> m2 s-1]. real :: Lat_fusion !< The latent heat of fusion [Q ~> J kg-1]. - real :: Gamma_T_3EQ !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation - real :: Gamma_S_3EQ !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation - !< This number should be specified by the user. + real :: Gamma_T_3EQ !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation [nondim] + real :: Gamma_S_3EQ !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation [nondim] + !< This number should be specified by the user. real :: col_mass_melt_threshold !< An ocean column mass below the iceshelf below which melting !! does not occur [R Z ~> kg m-2] logical :: mass_from_file !< Read the ice shelf mass from a file every dt @@ -168,6 +168,7 @@ module MOM_ice_shelf !! shelf dynamics will be initialized logical :: data_override_shelf_fluxes !< True if the ice shelf surface mass fluxes can be !! written using the data_override feature (only for MOSAIC grids) + logical :: data_override_melt !< True if overriding melt rate and heat flux logical :: override_shelf_movement !< If true, user code specifies the shelf movement !! instead of using the dynamic ice-shelf mode. logical :: isthermo !< True if the ice shelf can exchange heat and @@ -194,12 +195,14 @@ module MOM_ice_shelf real :: dTFr_dp !< Partial derivative of freezing temperature with !! pressure [C T2 R-1 L-2 ~> degC Pa-1] real :: Zeta_N !< The stability constant xi_N = 0.052 from Holland & Jenkins '99 - !! divided by the von Karman constant VK. Was 1/8. - real :: Vk !< Von Karman's constant - dimensionless - real :: Rc !< critical flux Richardson number. - logical :: buoy_flux_itt_bug !< If true, fixes buoyancy iteration bug - logical :: salt_flux_itt_bug !< If true, fixes salt iteration bug - real :: buoy_flux_itt_threshold !< Buoyancy iteration threshold for convergence + !! divided by the von Karman constant VK [nondim]. Was 1/8. + real :: Vk !< Von Karman's constant [nondim] + real :: Rc !< critical flux Richardson number [nondim] + logical :: ustar_from_vel_bugfix !< If true, fixes ustar from ocean velocity bug + logical :: buoy_flux_itt_bugfix !< If true, fixes buoyancy iteration bug + logical :: salt_flux_itt_bugfix !< If true, fixes salt iteration bug + real :: buoy_flux_tol !< Fractional buoyancy iteration tolerance for convergence [nondim] + real :: tideamp_scaling_factor !< Tideamp scaling factor in friction velocity calculation !>@{ Diagnostic handles integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & @@ -294,11 +297,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) !! This is computed as part of the ISOMIP diagnostics. real :: time_step !< Length of time over which these fluxes will be applied [T ~> s]. real :: Itime_step !< Inverse of the length of time over which these fluxes will be applied [T-1 ~> s-1] - real :: VK !< Von Karman's constant - dimensionless + real :: VK !< Von Karman's constant [nondim] real :: ZETA_N !< This is the stability constant xi_N = 0.052 from Holland & Jenkins '99 !! divided by the von Karman constant VK. Was 1/8. [nondim] - real :: RC !< critical flux Richardson number. - real :: I_ZETA_N !< The inverse of ZETA_N [nondim]. + real :: Rf_crit !< critical flux Richardson number [nondim] + real :: I_2Zeta_N !< Half the inverse of Zeta_N [nondim]. real :: I_LF !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1]. real :: I_VK !< The inverse of the Von Karman constant [nondim]. real :: PR, SC !< The Prandtl number and Schmidt number [nondim]. @@ -318,7 +321,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) real :: wB_flux !< The downward vertical flux of buoyancy just inside the ocean [Z2 T-3 ~> m2 s-3]. real :: dB_dS !< The derivative of buoyancy with salinity [Z T-2 S-1 ~> m s-2 ppt-1]. real :: dB_dT !< The derivative of buoyancy with temperature [Z T-2 C-1 ~> m s-2 degC-1]. - real :: I_n_star ! [nondim] + real :: I_n_star ! The inverse of the ratio of working boundary layer thickness + ! to the neutral thickness [nondim] real :: n_star_term ! A term in the expression for nstar [T3 Z-2 ~> s3 m-2] real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1] real :: dIns_dwB !< The partial derivative of I_n_star with wB_flux, in [T3 Z-2 ~> s3 m-2] @@ -327,34 +331,42 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) real :: dS_ustar ! The difference between the salinity at the ice-ocean interface and the ocean ! boundary layer salinity times the friction velocity [S Z T-1 ~> ppt m s-1] real :: ustar_h ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1] - real :: Gam_turb ! [nondim] + real :: Gam_turb ! A relative turbluent diffusivity [nondim] real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivities [nondim] real :: RhoCp ! A typical ocean density times the heat capacity of water [Q R C-1 ~> J m-3 degC-1] - real :: ln_neut + real :: ln_neut ! The log of the ratio of the neutral boundary layer thickness to the molecular + ! boundary layer thickness if it is greater than 1 or 0 otherwise [nondim] real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1] real :: Sb_min, Sb_max ! Minimum and maximum boundary salinities [S ~> ppt] real :: dS_min, dS_max ! Minimum and maximum salinity changes [S ~> ppt] ! Variables used in iterating for wB_flux. - real :: wB_flux_new, dDwB_dwB_in - real :: I_Gam_T, I_Gam_S - real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] + real :: wB_flux_next ! The next interation's guess for wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_new ! An updated value of wB_flux when Gam_turb is based on wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_max ! The upper bound on wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_min ! The lower bound on wB_flux [Z2 T-3 ~> m2 s-2] + real :: dDwB_dwB ! The slope of the change in wB_flux between iterations with wB_flux [nondim] + real :: DwB_max ! The change in wB_flux when it is wB_flux_max [Z2 T-3 ~> m2 s-2] + real :: DwB_min ! The change in wB_flux when it is wB_flux_min [Z2 T-3 ~> m2 s-2] + real :: I_Gam_T, I_Gam_S ! Terms that vary inversely with Gam_mol_T or Gam_mol_S and Gam_turb [nondim] + real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2]. real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2] - real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u- - real :: asv1, asv2 ! and v-points [L2 ~> m2]. + real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u-points [L2 ~> m2] + real :: asv1, asv2 ! Ocean areas covered by ice shelves at neighboring v-points [L2 ~> m2] real :: I_au, I_av ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2] real :: Irho0 ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1] logical :: Sb_min_set, Sb_max_set + logical :: root_found logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. logical :: coupled_GL ! If true, the grounding line position is determined based on ! coupled ice-ocean dynamics. - real, parameter :: c2_3 = 2.0/3.0 - character(len=160) :: mesg ! The text of an error message + real, parameter :: c2_3 = 2.0/3.0 ! Two thirds [nondim] + character(len=320) :: mesg ! The text of an error message integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, ied, jed, it1, it3 - real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [Z L2 ~> m3] - !for all ice sheets, Antarctica only, or Greenland only [Z L2 ~> m3] + real :: vaf0, vaf0_A, vaf0_G ! The previous volumes above floatation [Z L2 ~> m3] + ! for all ice sheets, Antarctica only, or Greenland only if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") @@ -394,8 +406,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ! useful parameters ZETA_N = CS%Zeta_N VK = CS%Vk - RC = CS%Rc - I_ZETA_N = 1.0 / ZETA_N + Rf_crit = CS%Rc + I_2Zeta_N = 0.5 / CS%Zeta_N I_LF = 1.0 / CS%Lat_fusion SC = CS%kv_molec/CS%kd_molec_salt PR = CS%kv_molec/CS%kd_molec_temp @@ -454,19 +466,26 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) tauy2 = (((asv1 * (sfc_state%tauy_shelf(i,J-1)**2)) + (asv2 * (sfc_state%tauy_shelf(i,J)**2)) ) * I_av) endif u2_av = (((asu1 * (sfc_state%u(I-1,j)**2)) + (asu2 * sfc_state%u(I,j)**2)) * I_au) - v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asu2 * sfc_state%v(i,J)**2)) * I_av) + if (CS%ustar_from_vel_bugfix) then + v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asv2 * sfc_state%v(i,J)**2)) * I_av) + else + v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asu2 * sfc_state%v(i,J)**2)) * I_av) + endif if ((taux2 + tauy2 > 0.0) .and. .not.CS%ustar_shelf_from_vel) then if (CS%ustar_max >= 0.0) then fluxes%ustar_shelf(i,j) = MIN(CS%ustar_max, MAX(CS%ustar_bg, US%L_to_Z * & - sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2))) + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%tideamp_scaling_factor * & + (CS%cdrag*CS%utide(i,j)**2)))) else fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & - sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2)) + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%tideamp_scaling_factor * & + (CS%cdrag*CS%utide(i,j)**2))) endif else ! Take care of the cases when taux_shelf is not set or not allocated. fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_TO_Z * & - sqrt(CS%cdrag*((u2_av + v2_av) + CS%utide(i,j)**2))) + sqrt(CS%cdrag*((u2_av + v2_av) + CS%tideamp_scaling_factor * & + (CS%utide(i,j)**2)))) endif else ! There is no shelf here. fluxes%ustar_shelf(i,j) = 0.0 @@ -502,11 +521,12 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) if (absf*sfc_state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = sfc_state%Hml(i,j) else ; hBL_neut = (VK*ustar_h) / absf ; endif hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%kv_molec)) + ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) + n_star_term = (ZETA_N * hBL_neut * VK) / (Rf_crit * ustar_h**3) ! Determine the mixed layer buoyancy flux, wB_flux. dB_dS = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dS(i) dB_dT = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dT(i) - ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) if (CS%find_salt_root) then ! Solve for the skin salinity using the linearized liquidus parameters and @@ -556,68 +576,152 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) dT_ustar = (ISS%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h dS_ustar = (Sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h - ! First, determine the buoyancy flux assuming no effects of stability - ! on the turbulence. Following H & J '99, this limit also applies - ! when the buoyancy flux is destabilizing. - - if (CS%const_gamma) then ! if using a constant gamma_T - ! note the different form, here I_Gam_T is NOT 1/Gam_T! + if (CS%const_gamma) then + ! If using a constant gamma_T, there are no effects of the buoyancy flux on the turbulence. I_Gam_T = CS%Gamma_T_3EQ I_Gam_S = CS%Gamma_S_3EQ - else - Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0)) + wT_flux = dT_ustar * CS%Gamma_T_3EQ + wB_flux = dB_dS * (dS_ustar * CS%Gamma_S_3EQ) + dB_dT * wT_flux + elseif (.not.CS%buoy_flux_itt_bugfix) then + ! Gamma_T and gamma_S are a function of the buoyancy flux, and there should have been + ! iteration to find the root where wB_flux is consistent with the values of gamma with + ! that flux, but it was omitted. + Gam_turb = I_VK * (ln_neut + (I_2Zeta_N - 1.0)) I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) - endif + wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T) - wT_flux = dT_ustar * I_Gam_T - wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux + if (wB_flux < 0.0) then ! The stabilising buoyancy flux reduces the turbulent fluxes. + I_n_star = sqrt(1.0 - n_star_term * wB_flux) + if (hBL_neut_h_molec > I_n_star**2) then + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + else ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + endif + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + endif + wT_flux = dT_ustar * I_Gam_T + else ! gamma_T and gamma_S are a function of the buoyancy flux with proper iteration. + ! Find the root where wB_flux is consistent with the values of gamma with that flux. + + ! First, determine the buoyancy flux assuming no effects of stability + ! on the turbulence. Following H & J '99, this limit also applies + ! when the buoyancy flux is destabilizing. + Gam_turb = I_VK * (ln_neut + (I_2Zeta_N - 1.0)) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T - if (wB_flux < 0.0) then - ! The buoyancy flux is stabilizing and will reduce the turbulent - ! fluxes, and iteration is required. - n_star_term = (ZETA_N * hBL_neut * VK) / (RC * ustar_h**3) - do it3 = 1,30 - ! n_star <= 1.0 is the ratio of working boundary layer thickness - ! to the neutral thickness. - ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL + if (wB_flux < 0.0) then + ! The buoyancy flux is stabilizing and will reduce the turbulent + ! fluxes, and iteration is required. + ! n_star <= 1.0 is the ratio of working boundary layer thickness + ! to the neutral thickness. I_n_star is its inverse. I_n_star = sqrt(1.0 - n_star_term * wB_flux) - dIns_dwB = 0.5 * n_star_term / I_n_star if (hBL_neut_h_molec > I_n_star**2) then - Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + & - (0.5*I_ZETA_N*I_n_star - 1.0)) - dG_dwB = I_VK * ( -2.0 / I_n_star + (0.5 * I_ZETA_N)) * dIns_dwB - else - ! The layer dominated by molecular viscosity is smaller than - ! the assumed boundary layer. This should be rare! - Gam_turb = I_VK * (0.5 * I_ZETA_N*I_n_star - 1.0) - dG_dwB = I_VK * (0.5 * I_ZETA_N) * dIns_dwB + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + else ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) endif - - if (CS%const_gamma) then ! if using a constant gamma_T - ! note the different form, here I_Gam_T is NOT 1/Gam_T! - I_Gam_T = CS%Gamma_T_3EQ - I_Gam_S = CS%Gamma_S_3EQ - else - I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) - I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + root_found = (abs(wB_flux_new - wB_flux) < CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) + ! Do not update the flux if its maagnitude would be increased by the otherwise + ! stabilizing buoyancy fluxes. This can happen when the buoyancy flux + ! is stabilizing when one of the heat or salt fluxes are destabilizing due + ! to their different molecular properties. + if (wB_flux_new <= wB_flux) root_found = .true. + + if (.not.root_found) then + wB_flux_max = 0.0 ; DwB_max = wB_flux + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + + if ((wB_flux_min*n_star_term < (1.0 - hBL_neut_h_molec)) .and. & + ((1.0 - hBL_neut_h_molec) < wB_flux_max*n_star_term)) then + ! The derivative of Gam_turb with wB_flux has a discontinuous change within the + ! bracketed range of values. Take this discontinous slope value for a first + ! guess, because Newton's method and the false position method may not converge + ! quickly when this discontinuity is between a guess and the solution. + wB_flux = (1.0 - hBL_neut_h_molec) / n_star_term + I_n_star = sqrt(hBL_neut_h_molec) + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + + if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) then + ! The root has been found to within the tolerance at the kink. This should be very rare. + root_found = .true. + elseif (wB_flux_new > wB_flux) then + ! The solution is in the limit where abs(wB_flux) is small and + ! Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + else + ! The solution is in the limt where abs(wB_flux) is large and + ! Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux + endif + endif endif - wT_flux = dT_ustar * I_Gam_T - wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux - - ! Find the root where wB_flux_new = wB_flux. - if (abs(wB_flux_new - wB_flux) < CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) exit + if (.not.root_found) then + ! Use the false position for the next guess. + wB_flux = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max)) + + do it3 = 1,30 + ! Iterate using Newton's method with bounds or the false position method to find the root. + + I_n_star = sqrt(1.0 - n_star_term * wB_flux) + dIns_dwB = -0.5 * n_star_term / I_n_star + if (hBL_neut_h_molec > I_n_star**2) then + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + dG_dwB = I_VK * (( -2.0 / I_n_star + I_2Zeta_N) * dIns_dwB) + else + ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + dG_dwB = I_VK * (I_2Zeta_N * dIns_dwB) + endif + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + + ! Test for convergence to within tolerance at the point where wB_flux_new = wB_flux. + if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) & + root_found = .true. + if (root_found) exit + + dDwB_dwB = -dG_dwB * ((dB_dS * dS_ustar) * I_Gam_S**2 + & + (dB_dT * dT_ustar) * I_Gam_T**2) - 1.0 + if ((dDwB_dwB >= 0.0) .or. & + ( wB_flux - wB_flux_new >= abs(dDwB_dwB)*(wB_flux_max - wB_flux)) .or. & + ( wB_flux - wB_flux_new <= abs(dDwB_dwB)*(wB_flux_min - wB_flux)) ) then + ! Use the False position method to determine the guess for the next iteration when + ! Newton's method would go out of bounds + wB_flux_next = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max)) + else + ! Use Newton's method for the next guess. + wB_flux_next = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB + endif + + ! Reset one of the bounds inward. + if (wB_flux_new - wB_flux > 0) then + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + else + wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux + endif + + ! Update wB_flux + wB_flux = wB_flux_next + enddo ! it3 + endif - dDwB_dwB_in = dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + & - dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0 - ! This is Newton's method without any bounds. Should bounds be needed? - wB_flux_new = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in - ! Update wB_flux - if (CS%buoy_flux_itt_bug) wB_flux = wB_flux_new - enddo !it3 - endif + endif ! End of test for first guess of wB_flux < 0. + wT_flux = dT_ustar * I_Gam_T + endif ! End of test for CS%const_gamma ISS%tflux_ocn(i,j) = RhoCp * wT_flux exch_vel_t(i,j) = ustar_h * I_Gam_T @@ -688,7 +792,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) Sbdry(i,j) = Sbdry_it endif ! Sb_min_set - if (.not.CS%salt_flux_itt_bug) Sbdry(i,j) = Sbdry_it + if (.not.CS%salt_flux_itt_bugfix) Sbdry(i,j) = Sbdry_it endif ! CS%find_salt_root @@ -1138,7 +1242,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step) type(ice_shelf_CS), pointer :: CS !< This module's control structure. type(surface), intent(inout) :: sfc_state !< Surface ocean state type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. - real, intent(in) :: time_step !< Time step over which fluxes are applied + real, intent(in) :: time_step !< Time step over which fluxes are applied [T ~> s] ! local variables real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. real :: frac_open !< The fractional area of the ocean that is not covered by the ice shelf [nondim]. @@ -1196,6 +1300,14 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step) enddo ; enddo endif + !CLAIRE + if (CS%data_override_melt) then + call data_override(G%Domain, 'water_flux', ISS%water_flux(is:ie,js:je), CS%Time, & + scale=US%kg_m2s_to_RZ_T) + call data_override(G%Domain, 'tflux_ocn', ISS%tflux_ocn(is:ie,js:je), CS%Time, & + scale=US%W_m2_to_QRZ_T) + endif + if (CS%debug) then call MOM_forcing_chksum("Before adding shelf fluxes", fluxes, G, CS%US, haloshift=0) endif @@ -1377,10 +1489,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, type(directories) :: dirs type(dyn_horgrid_type), pointer :: dG => NULL() type(dyn_horgrid_type), pointer :: dG_in => NULL() - real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic. + real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic + ! [T kg R-1 Z-1 m-2 s-1 ~> nondim] real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered ! to be floating when CONST_SEA_LEVEL = True [Z ~> m]. - real :: cdrag, drag_bg_vel + real :: cdrag ! The drag coefficient at the ice-ocean interface [nondim] + real :: drag_bg_vel ! A background velocity used in the quadratic drag [Z T-1 ~> m s-1] logical :: new_sim, save_IC !This include declares and sets the variable "version". # include "version_variable.h" @@ -1396,7 +1510,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, real :: utide ! A tidal velocity [L T-1 ~> m s-1] real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting ! does not occur [Z ~> m] - real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data + real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for ice shelf input data [L T-1 ~> m s-1] type(surface), pointer :: sfc_state => NULL() type(vardesc) :: u_desc, v_desc @@ -1540,6 +1654,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, if (CS%GL_regularize) CS%GL_couple = .false. if (CS%solo_ice_sheet) CS%GL_couple = .false. endif + call get_param(param_file, mdl, "DATA_OVERRIDE_MELT", & + CS%data_override_melt, & + "If true, the data override feature is used for melt rate "//& + "and heat flux from ice shelves.", default=.false.) call get_param(param_file, mdl, "SHELF_THERMO", CS%isthermo, & "If true, use a thermodynamically interactive ice shelf.", & @@ -1686,11 +1804,14 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call get_param(param_file, mdl, "ICE_SHELF_RC", CS%Rc, & "Critical flux Richardson number for ice melt ", & units="nondim", default=0.20) - call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUG", CS%buoy_flux_itt_bug, & - "Bug fix of buoyancy iteration", default=.true.) - call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUG", CS%salt_flux_itt_bug, & - "Bug fix of salt iteration", default=.true.) - call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_itt_threshold, & + call get_param(param_file, mdl, "ICE_SHELF_USTAR_FROM_VEL_BUGFIX", CS%ustar_from_vel_bugfix, & + "Bug fix for ice-area weighting of squared ocean velocities "//& + "used to calculate friction velocity under ice shelves", default=.false.) + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX", CS%buoy_flux_itt_bugfix, & + "Bug fix of buoyancy iteration", default=.true., old_name="ICE_SHELF_BUOYANCY_FLUX_ITT_BUG") + call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUGFIX", CS%salt_flux_itt_bugfix, & + "Bug fix of salt iteration", default=.true., old_name="ICE_SHELF_SALT_FLUX_ITT_BUG") + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_tol, & "Convergence criterion of Newton's method for ice shelf "//& "buoyancy iteration.", units="nondim", default=1.0e-4) @@ -1709,6 +1830,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 + call get_param(param_file, mdl, "ICE_SHELF_TIDEAMP_SCALING_FACTOR", CS%tideamp_scaling_factor, & + "Scale TIDEAMP_FILE or UTIDE by number in melt parameterisation "//& + "friction velocity calculation.", units = "nondim", default=1.0) if (read_TIDEAMP) then call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & "The path to the file containing the spatially varying tidal amplitudes.", & @@ -1903,7 +2027,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, elseif (.not.new_sim) then ! This line calls a subroutine that reads the initial conditions from a restart file. call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.") - call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp) + ! call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp) + call restore_state('Shelf.res.nc', dirs%restart_input_dir, Time, G, CS%restart_CSp) endif ! .not. new_sim @@ -2356,6 +2481,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) end select end subroutine initialize_shelf_mass + !> This subroutine applies net accumulation/ablation at the top surface to the dynamic ice shelf. !! acc_rate[m-s]=surf_mass_flux/density_ice is ablation/accumulation rate !! positive for accumulation negative for ablation @@ -2372,14 +2498,13 @@ subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time ! locals integer :: i, j - real ::I_rho_ice + real :: I_rho_ice ! The specific volume of ice [R-1 ~> m3 kg-1] I_rho_ice = 1.0 / CS%density_ice !update time ! CS%Time = Time - ! CS%time_step = time_step ! update surface mass flux rate ! if (CS%surf_mass_flux_from_file) call update_surf_mass_flux(G, US, CS, ISS, Time) @@ -2459,13 +2584,14 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) end subroutine update_shelf_mass !> Save the ice shelf restart file -subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes) +subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes, data_override_melt) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nondim]. - real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf ! kg m-2] + real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2] logical, optional :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using !! the data_override capability (only for MOSAIC grids) + logical, optional :: data_override_melt !< If true, basal melt overridden integer :: i, j @@ -2488,6 +2614,11 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_ if (CS%active_shelf_dynamics) data_override_shelf_fluxes = CS%data_override_shelf_fluxes endif + if (present(data_override_melt)) then + data_override_melt = .false. + if (CS%data_override_melt) data_override_melt = CS%data_override_melt + endif + end subroutine ice_shelf_query !> Save the ice shelf restart file diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 9f50a77881..be1a08919b 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -2540,7 +2540,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif if (CS%max_surface_slope>0) then - scale = min(CS%max_surface_slope/sqrt((sx**2)+(sy**2)),1.0) + scale = CS%max_surface_slope / max( sqrt((sx**2) + (sy**2)), CS%max_surface_slope ) sx = scale*sx; sy = scale*sy endif diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index b6d4dfa489..d285d08b94 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -107,7 +107,7 @@ module MOM_diabatic_aux !! This subroutine warms any water that is colder than the (currently !! surface) freezing point up to the freezing point and accumulates !! the required heat (in [Q R Z ~> J m-2]) in tv%frazil. -subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) +subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo, frac_shelf_h) 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)), & @@ -120,6 +120,8 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]. integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage ! Local variables real, dimension(SZI_(G)) :: & @@ -188,7 +190,39 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif endif ; enddo endif + if (PRESENT(frac_shelf_h)) then + do k=nz,1,-1 + T_fr_set = .false. + do i=is,ie + if (((G%mask2dT(i,j) > 0.0) .and. (.not. (frac_shelf_h(i,j) > 0.0))) .and. & + ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then + if (.not.T_fr_set) then + call calculate_TFreeze(tv%S(i:ie,j,k), pressure(i:ie,k), T_freeze(i:ie), & + tv%eqn_of_state) + T_fr_set = .true. + endif + + hc = (tv%C_p*GV%H_to_RZ) * h(i,j,k) + if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then + ! Very thin layers should not be cooled by the frazil flux. + if (tv%T(i,j,k) < T_freeze(i)) then + fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) + tv%T(i,j,k) = T_freeze(i) + endif + elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then + if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then + tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc + fraz_col(i) = 0.0 + else + fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) + tv%T(i,j,k) = T_freeze(i) + endif + endif + endif + enddo + enddo + else do k=nz,1,-1 T_fr_set = .false. do i=is,ie @@ -219,6 +253,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif enddo enddo + endif do i=is,ie tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i) enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 95c1d3d265..63bcbc4033 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -183,6 +183,8 @@ module MOM_diabatic_driver logical :: Use_KdWork_diag = .false. !< Logical flag to indicate if any Kd_work diagnostics are on. logical :: Use_N2_diag = .false. !< Logical flag to indicate if any N2 diagnostics are on. + logical :: frazil_not_under_iceshelves !< If true, turn off frazil under ice shelves + !>@{ Diagnostic IDs integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 @@ -374,7 +376,12 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + if (CS%frazil_not_under_iceshelves) then + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff, & + frac_shelf_h=fluxes%frac_shelf_h) + else + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + endif else call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, halo=CS%halo_TS_diff) endif @@ -434,7 +441,12 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full) + if (CS%frazil_not_under_iceshelves) then + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff, & + frac_shelf_h=fluxes%frac_shelf_h) + else + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + endif else call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp) endif @@ -3686,6 +3698,10 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di endif endif + call get_param(param_file, mdl, "FRAZIL_NOT_UNDER_ICESHELF", CS%frazil_not_under_iceshelves, & + "If true, do not use frazil scheme underneath ice shelves "//& + "defined by frac_shelf_h greater than 0.", default=.false.) + ! diagnostics for tendencies of temp and heat due to frazil CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, & long_name='Cell Thickness', standard_name='cell_thickness', & diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 new file mode 100644 index 0000000000..e092aa7d6d --- /dev/null +++ b/src/tracer/MOM_generic_tracer.F90 @@ -0,0 +1,946 @@ +!> Drives the generic version of tracers TOPAZ and CFC and other GFDL BGC components +module MOM_generic_tracer + +! This file is part of MOM6. See LICENSE.md for the license. + +#include + +! The following macro is usually defined in but since MOM6 should not directly +! include files from FMS we replicate the macro lines here: +#ifdef NO_F2000 +#define _ALLOCATED associated +#else +#define _ALLOCATED allocated +#endif + + ! ### These imports should not reach into FMS directly ### + use field_manager_mod, only: fm_string_len + + use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list + use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag + use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set + use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer + use generic_tracer, only: generic_tracer_update_from_bottom,generic_tracer_vertdiff_G + use generic_tracer, only: generic_tracer_coupler_accumulate, generic_tracer_update_from_coupler + + use g_tracer_utils, only: g_tracer_get_name,g_tracer_set_values,g_tracer_set_common,g_tracer_get_common + use g_tracer_utils, only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init + use g_tracer_utils, only: g_tracer_send_diag,g_tracer_get_values + use g_tracer_utils, only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag + use g_tracer_utils, only: g_tracer_get_obc_segment_props + + use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS + use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here + use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr + use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end + use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe + use MOM_file_parser, only : get_param, log_param, log_version, param_file_type + use MOM_forcing_type, only : forcing, optics_type + use MOM_grid, only : ocean_grid_type + use MOM_hor_index, only : hor_index_type + use MOM_interface_heights, only : thickness_to_dz + use MOM_io, only : file_exists, MOM_read_data, slasher + use MOM_open_boundary, only : ocean_OBC_type + use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments + use MOM_open_boundary, only : set_obgc_segments_props + use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS + use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP, array_global_min_max + use MOM_sponge, only : set_up_sponge_field, sponge_CS + use MOM_time_manager, only : time_type, set_time + use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut + use MOM_tracer_registry, only : register_tracer, tracer_registry_type + use MOM_tracer_Z_init, only : tracer_Z_init + use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z + use MOM_unit_scaling, only : unit_scale_type + use MOM_variables, only : surface, thermo_var_ptrs + use MOM_verticalGrid, only : verticalGrid_type + + + implicit none ; private + + !> A state hidden in module data that is very much not allowed in MOM6 + ! ### This needs to be fixed + logical :: g_registered = .false. + + public register_MOM_generic_tracer, initialize_MOM_generic_tracer + public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state + public end_MOM_generic_tracer, MOM_generic_tracer_get + public MOM_generic_tracer_stock + public MOM_generic_flux_init + public MOM_generic_tracer_min_max + public MOM_generic_tracer_fluxes_accumulate + public register_MOM_generic_tracer_segments + + !> Control structure for generic tracers + type, public :: MOM_generic_tracer_CS ; private + character(len = 200) :: IC_file !< The file in which the generic tracer initial values can + !! be found, or an empty string for internal initialization. + logical :: Z_IC_file !< If true, the generic_tracer IC_file is in Z-space. The default is false. + real :: tracer_IC_val = 0.0 !< The initial value assigned to tracers, in + !! concentration units [conc] + real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out, in + !! concentration units [conc] + logical :: tracers_may_reinit !< If true, tracers may go through the + !! initialization code if they are not found in the restart files. + + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Restart control structure + type(ocean_OBC_type), pointer :: OBC => NULL() ! Pointer to the first element of the linked list of generic tracers. + type(g_tracer_type), pointer :: g_tracer_list => NULL() + + end type MOM_generic_tracer_CS + +contains + + !> Initializes the generic tracer packages and adds their tracers to the list + !! Adds the tracers in the list of generic tracers to the set of MOM tracers (i.e., MOM-register them) + !! Register these tracers for restart + function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< Horizontal index ranges + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct + + ! Local variables + logical :: register_MOM_generic_tracer + logical :: obc_has + ! This include declares and sets the variable "version". +# include "version_variable.h" + + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' + character(len=200) :: inputdir ! The directory where NetCDF input files are. + ! These can be overridden later in via the field manager? + + integer :: ntau, axes(3) + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + character(len=fm_string_len) :: g_tracer_name, longname,units + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(HI),SZJ_(HI)) :: grid_kmt ! A 2-d array of nk + + register_MOM_generic_tracer = .false. + if (associated(CS)) then + call MOM_error(FATAL, "register_MOM_generic_tracer called with an "// & + "associated control structure.") + endif + allocate(CS) + + + !Register all the generic tracers used and create the list of them. + !This can be called by ALL PE's. No array fields allocated. + if (.not. g_registered) then + call generic_tracer_register() + g_registered = .true. + endif + + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, sub_name, version, "") + call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, & + "The file in which the generic tracer initial values can "//& + "be found, or an empty string for internal initialization.", & + default=" ") + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + ! Add the directory if CS%IC_file is not already a complete path. + call get_param(param_file, sub_name, "INPUTDIR", inputdir, default=".") + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, sub_name, "INPUTDIR/GENERIC_TRACER_IC_FILE", CS%IC_file) + endif + call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE_IS_Z", CS%Z_IC_file, & + "If true, GENERIC_TRACER_IC_FILE is in depth space, not "//& + "layer space.",default=.false.) + call get_param(param_file, sub_name, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + + CS%restart_CSp => restart_CS + + ntau=1 ! MOM needs the fields at only one time step + + + ! At this point G%mask2dT and CS%diag%axesTL are not allocated. + ! postpone diag_registeration to initialize_MOM_generic_tracer + + !Fields cannot be diag registered as they are allocated and have to registered later. + grid_tmask(:,:,:) = 0.0 + grid_kmt(:,:) = 0 + axes(:) = -1 + + ! + ! Initialize all generic tracers + ! + call generic_tracer_init(HI%isc,HI%iec,HI%jsc,HI%jec,HI%isd,HI%ied,HI%jsd,HI%jed,& + GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0)) + + + ! + ! MOM-register the generic tracers + ! + + !Get the tracer list + call generic_tracer_get_list(CS%g_tracer_list) + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + ! For each tracer name get its T_prog index and get its fields + + g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,g_tracer_name) + + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) + call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname) + call g_tracer_get_values(g_tracer,g_tracer_name,'units',units ) + + !!nnz: MOM field is 3D. Does this affect performance? Need it be override field? + tr_ptr => tr_field(:,:,:,1) + ! Register prognostic tracer for horizontal advection, diffusion, and restarts. + if (g_tracer_is_prog(g_tracer)) then + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + name=g_tracer_name, longname=longname, units=units, & + registry_diags=.false., & !### CHANGE TO TRUE? + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + else + call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & + restart_CS, longname=longname, units=units) + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + + register_MOM_generic_tracer = .true. + end function register_MOM_generic_tracer + + !> Register OBC segments for generic tracers + subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file) + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + + ! Local variables + logical :: obc_has + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments' + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + character(len=fm_string_len) :: g_tracer_name + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] + + if (.NOT. associated(OBC)) return + !Get the tracer list + call generic_tracer_get_list(CS%g_tracer_list) + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + + g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,g_tracer_name) + if (g_tracer_is_prog(g_tracer)) then + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,& + obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + if (obc_has) then + call set_obgc_segments_props(OBC,g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + call register_obgc_segments(GV, OBC, tr_Reg, param_file, g_tracer_name) + endif + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + + end subroutine register_MOM_generic_tracer_segments + + !> Initialize phase II: Initialize required variables for generic tracers + !! There are some steps of initialization that cannot be done in register_MOM_generic_tracer + !! This is the place and time to do them: + !! Set the grid mask and initial time for all generic tracers. + !! Diag_register them. + !! Z_diag_register them. + !! + !! This subroutine initializes the NTR tracer fields in tr(:,:,:,:) + !! and it sets up the tracer output. + subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, & + CS, sponge_CSp, ALE_sponge_CSp) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic + !! variables + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(sponge_CS), pointer :: sponge_CSp !< Pointer to the control structure for the sponges. + type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the + !! ALE sponges. + + character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer' + logical :: OK,obc_has + integer :: i, j, k, isc, iec, jsc, jec, nk + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + character(len=fm_string_len) :: g_tracer_name + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer vertical extent [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(G),SZJ_(G)) :: grid_kmt ! A 2-d array of nk + + !! 2010/02/04 Add code to re-initialize Generic Tracers if needed during a model simulation + !! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped. + !! Ideally, the generic tracer IC file should have the tracers on Z levels. + + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke + + CS%diag=>diag + !Get the tracer list + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + !For each tracer name get its fields + g_tracer=>CS%g_tracer_list + + call thickness_to_dz(h, tv, dz, G, GV, US) + + do + if (INDEX(CS%IC_file, '_NULL_') /= 0) then + call MOM_error(WARNING, "The name of the IC_file "//trim(CS%IC_file)//& + " indicates no MOM initialization was asked for the generic tracers."//& + "Bypassing the MOM initialization of ALL generic tracers!") + exit + endif + call g_tracer_get_alias(g_tracer,g_tracer_name) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) + tr_ptr => tr_field(:,:,:,1) + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then + + if (g_tracer%requires_src_info ) then + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initializing generic tracer "//trim(g_tracer_name)//& + " using MOM_initialize_tracer_from_Z ") + + call MOM_initialize_tracer_from_Z(dz, tr_ptr, G, GV, US, param_file, & + src_file=g_tracer%src_file, src_var_nam=g_tracer%src_var_name, & + src_var_unit_conversion=g_tracer%src_var_unit_conversion, & + src_var_record=g_tracer%src_var_record, src_var_gridspec=g_tracer%src_var_gridspec, & + h_in_Z_units=.true.) + + !Check/apply the bounds for each g_tracer + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min + !Jasmin does not want to apply the maximum for now + !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max + endif + enddo ; enddo ; enddo + + !jgj: Reset CASED to 0 below K=1 + if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then + do k=2,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + tr_ptr(i,j,k) = 0.0 + endif + enddo ; enddo ; enddo + endif + elseif(.not. g_tracer%requires_restart) then + !Do nothing for this tracer, it is initialized by the tracer package + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "skip initialization of generic tracer "//trim(g_tracer_name)) + else !Do it old way if the tracer is not registered to start from a specific source file. + !This path should be deprecated if all generic tracers are required to start from specified sources. + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, & + "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//& + "Unable to read "//trim(g_tracer_name)//" from "//& + trim(CS%IC_file)//".") + endif + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initialized generic tracer "//trim(g_tracer_name)//& + " using Generic Tracer File on Z: "//CS%IC_file) + else + ! native grid + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) + call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain) + endif + else + call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//& + "check Generic Tracer IC filename "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) + endif + + endif + + call set_initialized(tr_ptr, g_tracer_name, CS%restart_CSp) + endif + + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ) + if(obc_has .and. g_tracer_is_prog(g_tracer) .and. .not.restart) & + call fill_obgc_segments(G, GV, OBC, tr_ptr, g_tracer_name) + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + enddo + !! end section to re-initialize generic tracers + + + !Now we can reset the grid mask, axes and time to their true values + !Note that grid_tmask must be set correctly on the data domain boundary + !so that coast mask can be deduced from it. + grid_tmask(:,:,:) = 0.0 + grid_kmt(:,:) = 0 + do j = G%jsd, G%jed ; do i = G%isd, G%ied + if (G%mask2dT(i,j) > 0.0) then + grid_tmask(i,j,:) = 1.0 + grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer. + endif + enddo ; enddo + call g_tracer_set_common(G%isc,G%iec,G%jsc,G%jec,G%isd,G%ied,G%jsd,G%jed,& + GV%ke,1,CS%diag%axesTL%handles,grid_tmask,grid_kmt,day) + + ! Register generic tracer modules diagnostics + +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + call generic_tracer_register_diag() +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + + end subroutine initialize_MOM_generic_tracer + + !> Column physics for generic tracers. + !! Get the coupler values for generic tracers that exchange with atmosphere + !! Update generic tracer concentration fields from sources and sinks. + !! Vertically diffuse generic tracer concentration fields. + !! Update generic tracers from bottom and their bottom reservoir. + !! + !! This subroutine applies diapycnal diffusion and any other column + !! tracer physics or chemistry to the tracers from this file. + !! CFCs are relatively simple, as they are passive tracers. with only a surface + !! flux as a source. + subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, US, CS, tv, optics, & + evap_CFL_limit, minimum_forcing_depth) + 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(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: ea !< The amount of fluid entrained from the layer + !! above during this call [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: eb !< The amount of fluid entrained from the layer + !! below during this call [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic + !! and tracer forcing fields. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml !< Mixed layer depth [Z ~> m] + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(optics_type), intent(in) :: optics !< The structure containing optical properties. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + ! Stored previously in diabatic CS. + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes + !! can be applied [H ~> m or kg m-2] + ! Stored previously in diabatic CS. + ! The arguments to this subroutine are redundant in that + ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + + ! Local variables + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics' + + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + character(len=fm_string_len) :: g_tracer_name + real, dimension(:,:), pointer :: stf_array ! The surface flux of the tracer [conc kg m-2 s-1] + real, dimension(:,:), pointer :: trunoff_array ! The tracer concentration in the river runoff [conc] + real, dimension(:,:), pointer :: runoff_tracer_flux_array ! The runoff tracer flux [conc kg m-2 s-1] + + real :: surface_field(SZI_(G),SZJ_(G)) ! The surface value of some field, here only used for salinity [S ~> ppt] + real :: dz_ml(SZI_(G),SZJ_(G)) ! The mixed layer depth in the MKS units used for generic tracers [m] + real :: sosga ! The global mean surface salinity [ppt] + + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! A work array of thicknesses [H ~> m or kg m-2] + integer :: i, j, k, isc, iec, jsc, jec, nk + + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke + + !Get the tracer list + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,& + trim(sub_name)//": No tracer in the list.") + +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + + ! + !Extract the tracer surface fields from coupler and update tracer fields from sources + ! + !call generic_tracer_coupler_get(fluxes%tr_fluxes) + !Niki: This is moved out to ocean_model_MOM.F90 because if dt_therm>dt_cpld we need to average + ! the fluxes without coming into this subroutine. + ! MOM5 has to modified to conform. + + ! + !Call the generic_tracer's update_from_coupler routine (convert salt_flux_added to g/m^2/sec) + ! + call generic_tracer_update_from_coupler(G%isd, G%jsd, 1000*(US%RZ_T_to_kg_m2s*fluxes%salt_flux_added)) + + ! + !Add contribution of river to surface flux + ! + g_tracer=>CS%g_tracer_list + do + if (_ALLOCATED(g_tracer%trunoff) .and. (.NOT. g_tracer%runoff_added_to_stf)) then + call g_tracer_get_alias(g_tracer,g_tracer_name) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array) + !nnz: Why is fluxes%river = 0? + runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * & + US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:) + stf_array = stf_array + runoff_tracer_flux_array + g_tracer%runoff_added_to_stf = .true. + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer => g_tracer_next + + enddo + + ! + !Prepare input arrays for source update + ! + + rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom_H + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k) + enddo ; enddo ; enddo + + dzt(:,:,:) = 1.0 + call thickness_to_dz(h_old, tv, dzt, G, GV, US) + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + dzt(i,j,k) = US%Z_to_m * dzt(i,j,k) + enddo ; enddo ; enddo + dz_ml(:,:) = 0.0 + do j=jsc,jec ; do i=isc,iec + surface_field(i,j) = tv%S(i,j,1) + dz_ml(i,j) = US%Z_to_m * Hml(i,j) + enddo ; enddo + sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt) + + ! + !Calculate tendencies (i.e., field changes at dt) from the sources / sinks + ! + if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. & + (G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. & + (US%C_to_degC == 1.0) .and. (US%S_to_ppt == 1.0)) then + ! Avoid unnecessary copies when no unit conversion is needed. + call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%areaT, get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & + internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga) + else + ! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True, + ! so we have to check and only do the scaling if it is associated. + if(associated(tv%internal_heat)) then + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & + internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), & + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga) + else + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga) + endif + endif + + ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes + ! usually in ALE mode + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + g_tracer=>CS%g_tracer_list + do + if (g_tracer_is_prog(g_tracer)) then + do k=1,nk ;do j=jsc,jec ; do i=isc,iec + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, & + fluxes, h_work, evap_CFL_limit, minimum_forcing_depth) + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + enddo + endif + + ! + !Update Tr(n)%field from explicit vertical diffusion + ! + ! Use a tridiagonal solver to determine the concentrations after the + ! surface source is applied and diapycnal advection and diffusion occurs. + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + ! Last arg is tau which is always 1 for MOM6 + call generic_tracer_vertdiff_G(h_work, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1) + else + ! Last arg is tau which is always 1 for MOM6 + call generic_tracer_vertdiff_G(h_old, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1) + endif + + ! Update bottom fields after vertical processes + + ! Second arg is tau which is always 1 for MOM6 + call generic_tracer_update_from_bottom(US%T_to_s*dt, 1, get_diag_time_end(CS%diag)) + + !Output diagnostics via diag_manager for all generic tracers and their fluxes + call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + + end subroutine MOM_generic_tracer_column_physics + + !> This subroutine calculates mass-weighted integral on the PE either + !! of all available tracer concentrations, or of a tracer that is + !! being requested specifically, returning the number of stocks it has + !! calculated. If the stock_index is present, only the stock corresponding + !! to that coded index is returned. + function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index) + 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(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(EFP_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each + !! tracer, in kg times concentration units [kg conc] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< The coded index of a specific stock + !! being sought. + integer :: MOM_generic_tracer_stock !< Return value, the + !! number of stocks calculated here. + + ! Local variables + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock' + + integer :: m + + MOM_generic_tracer_stock = 0 + if (.not.associated(CS)) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + if (.NOT. associated(CS%g_tracer_list)) return ! No stocks. + + m=1 ; g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,names(m)) + call g_tracer_get_values(g_tracer,names(m),'units',units(m)) + units(m) = trim(units(m))//" kg" + call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field) + + tr_ptr => tr_field(:,:,:,1) + stocks(m) = global_mass_int_EFP(h, G, GV, tr_ptr, on_PE_only=.true.) + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + m = m+1 + enddo + + MOM_generic_tracer_stock = m + + end function MOM_generic_tracer_stock + + !> This subroutine finds the global min and max of either of all available + !! tracer concentrations, or of a tracer that is being requested specifically, + !! returning the number of tracers it has evaluated. + !! It also optionally returns the locations of the extrema. + function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, names, units, & + xgmin, ygmin, zgmin, xgmax, ygmax, zgmax) + integer, intent(in) :: ind_start !< The index of the tracer to start with + logical, dimension(:), intent(out) :: got_minmax !< Indicates whether the global min and + !! max are found for each tracer + real, dimension(:), intent(out) :: gmin !< Global minimum of each tracer [conc] + real, dimension(:), intent(out) :: gmax !< Global maximum of each tracer [conc] + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + real, dimension(:), optional, intent(out) :: xgmin !< The x-position of the global minimum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, dimension(:), optional, intent(out) :: ygmin !< The y-position of the global minimum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, dimension(:), optional, intent(out) :: zgmin !< The z-position of the global minimum [layer] + real, dimension(:), optional, intent(out) :: xgmax !< The x-position of the global maximum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, dimension(:), optional, intent(out) :: ygmax !< The y-position of the global maximum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, dimension(:), optional, intent(out) :: zgmax !< The z-position of the global maximum [layer] + integer :: MOM_generic_tracer_min_max !< Return value, the + !! number of tracers done here. + + ! Local variables + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + real, dimension(:,:,:,:), pointer :: tr_field ! The tracer array whose extrema are being sought [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! The tracer array whose extrema are being sought [conc] + real :: x_min ! The x-position of the global minimum in the units of G%geoLonT, often [degrees_E] or [km] or [m] + real :: y_min ! The y-position of the global minimum in the units of G%geoLatT, often [degrees_N] or [km] or [m] + real :: z_min ! The z-position of the global minimum [layer] + real :: x_max ! The x-position of the global maximum in the units of G%geoLonT, often [degrees_E] or [km] or [m] + real :: y_max ! The y-position of the global maximum in the units of G%geoLatT, often [degrees_N] or [km] or [m] + real :: z_max ! The z-position of the global maximum [layer] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max' + + logical :: find_location + integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau + integer :: k, is, ie, js, je, m + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + MOM_generic_tracer_min_max = 0 + if (.not.associated(CS)) return + + if (.NOT. associated(CS%g_tracer_list)) return ! No stocks. + + call g_tracer_get_common(isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau) + find_location = present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. & + present(xgmax) .or. present(ygmax) .or. present(zgmax) + + m=ind_start ; g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,names(m)) + call g_tracer_get_values(g_tracer,names(m),'units',units(m)) + call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field) + + gmin(m) = -1.0 + gmax(m) = -1.0 + + tr_ptr => tr_field(:,:,:,1) + + if (find_location) then + call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m), & + x_min, y_min, z_min, x_max, y_max, z_max) + if (present(xgmin)) xgmin(m) = x_min + if (present(ygmin)) ygmin(m) = y_min + if (present(zgmin)) zgmin(m) = z_min + if (present(xgmax)) xgmax(m) = x_max + if (present(ygmax)) ygmax(m) = y_max + if (present(zgmax)) zgmax(m) = z_max + else + call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m)) + endif + + got_minmax(m) = .true. + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + m = m+1 + enddo + + MOM_generic_tracer_min_max = m + + end function MOM_generic_tracer_min_max + + !> This subroutine calculates the surface state and sets coupler values for + !! those generic tracers that have flux exchange with atmosphere. + !! + !! This subroutine sets up the fields that the coupler needs to calculate the + !! CFC fluxes between the ocean and atmosphere. + subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + ! Local variables + real :: sosga ! The global mean surface salinity [ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV),1) :: rho0 ! An unused array of densities [kg m-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] + + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state' + + !Set coupler values + !nnz: fake rho0 + rho0(:,:,:,:) = 1.0 + + dzt(:,:,:) = GV%H_to_m * h(:,:,:) + + sosga = global_area_mean(sfc_state%SSS, G, unscale=G%US%S_to_ppt) + + if ((G%US%C_to_degC == 1.0) .and. (G%US%S_to_ppt == 1.0)) then + call generic_tracer_coupler_set(sfc_state%tr_fields, & + ST=sfc_state%SST, SS=sfc_state%SSS, & + rho=rho0, & !nnz: required for MOM5 and previous versions. + ilb=G%isd, jlb=G%jsd, & + dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars + tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag)) + else + call generic_tracer_coupler_set(sfc_state%tr_fields, & + ST=G%US%C_to_degC*sfc_state%SST, SS=G%US%S_to_ppt*sfc_state%SSS, & + rho=rho0, & !nnz: required for MOM5 and previous versions. + ilb=G%isd, jlb=G%jsd, & + dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars + tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag)) + endif + + !Output diagnostics via diag_manager for all tracers in this module +! if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& +! "No tracer in the list.") +! call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) + !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld + ! hence if dt_therm > dt_cpld we get output (and contribution to the mean) at times that tracers + ! had not been updated. + ! Moving this to the end of column physics subroutine fixes this issue. + + end subroutine MOM_generic_tracer_surface_state + +!ALL PE subroutine on Ocean! Due to otpm design the fluxes should be initialized like this on ALL PE's! + subroutine MOM_generic_flux_init(verbosity) + integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity. + + character(len=128), parameter :: sub_name = 'MOM_generic_flux_init' + type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next + + if (.not. g_registered) then + call generic_tracer_register() + g_registered = .true. + endif + + call generic_tracer_get_list(g_tracer_list) + if (.NOT. associated(g_tracer_list)) then + call MOM_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") + return + endif + + g_tracer=>g_tracer_list + do + + call g_tracer_flux_init(g_tracer, verbosity=verbosity) + + ! traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + + end subroutine MOM_generic_flux_init + + subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight) + type(forcing), intent(in) :: flux_tmp !< A structure containing pointers to + !! thermodynamic and tracer forcing fields. + real, intent(in) :: weight !< A weight for accumulating this flux [nondim] + + call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight) + + end subroutine MOM_generic_tracer_fluxes_accumulate + + !> Copy the requested tracer into an array. + subroutine MOM_generic_tracer_get(name,member,array, CS) + character(len=*), intent(in) :: name !< Name of requested tracer. + character(len=*), intent(in) :: member !< The tracer element to return. + real, dimension(:,:,:), intent(out) :: array !< Array filled by this routine, in arbitrary units [A] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + ! Local variables + real, dimension(:,:,:), pointer :: array_ptr ! The tracer in the generic tracer structures, in + ! arbitrary units [A] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' + + call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr) + array(:,:,:) = array_ptr(:,:,:) + + end subroutine MOM_generic_tracer_get + + !> This subroutine deallocates the memory owned by this module. + subroutine end_MOM_generic_tracer(CS) + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + call generic_tracer_end() + + if (associated(CS)) then + deallocate(CS) + endif + end subroutine end_MOM_generic_tracer + +!---------------------------------------------------------------- +! Niki Zadeh +! +! +! William Cooke +! +! +! +! This module drives the generic version of tracers TOPAZ and CFC +! +!---------------------------------------------------------------- + +end module MOM_generic_tracer