From 857fcda6d89c08a5a734950321ee5a17d1684535 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Fri, 3 Nov 2023 10:23:36 +1100 Subject: [PATCH] allow generic tracers in NUOPC cap update build for generic_tracers initialise fms initialise gas exchange structures spawn Ice_ocean_boundary%fluxes register generic tracer coupler restarts write generic tracer flux restarts to nuopc-configured additional_restart_dir add function for getting CMEPS field name from FMS flux name copy atmos_ocean_fluxes_calc.F90 and atmos_ocean_dep_fluxes_calc.F90 from FMScoupler (commit: eeadda8dc74f605e9cb3b03be8a43ab2ce698496) update atmos_ocean_fluxes_calc.F90 (operate on 2D inputs, rather than 1D, add calculation for air_sea_deposition from 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, rename gas_fields_ice to gas_fields_ocn) move atmos_ocean_fluxes_calc routine into MOM_cap_gtracer_flux module add param array into Ice_ocean_boundary%fluxes after spawning add ability to override atm fields and IOB fluxes from data_table pass optional atm_fields to mom_import as 2D BC type do atm_fields override and flux calcs from mom_cap allow diagnostics for tracer fluxes and related fields spawn 2D coupler_type for atmos fields in InitializeAdvertise and set in ESMF internal state manually override IOB%fluxes since FMS coupler_type routine doesnt set override flag set u10 and psurf on atm_fields even when pcair cmeps field name is unknown to allow flux calculation with overridden pcair fields --- .gitignore | 5 + .gitmodules | 3 + MOM6/CMakeLists.txt | 58 +- MOM6/GFDL_generic_tracers | 1 + MOM6/extra_sources/mom_cap_gtracer_flux.F90 | 674 +++++++++++++++++++ MOM6/patches/MOM_coupler_types.F90.patch | 28 + MOM6/patches/MOM_couplertype_infra.F90.patch | 28 + MOM6/patches/mom_cap.F90.patch | 313 ++++++++- MOM6/patches/mom_cap_methods.F90.patch | 96 +++ MOM6/patches/mom_ocean_model_nuopc.F90.patch | 19 +- 10 files changed, 1205 insertions(+), 20 deletions(-) create mode 160000 MOM6/GFDL_generic_tracers create mode 100644 MOM6/extra_sources/mom_cap_gtracer_flux.F90 create mode 100644 MOM6/patches/MOM_coupler_types.F90.patch create mode 100644 MOM6/patches/MOM_couplertype_infra.F90.patch create mode 100644 MOM6/patches/mom_cap_methods.F90.patch diff --git a/.gitignore b/.gitignore index 1c66c69..89e8262 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,8 @@ build/* bin/* +Debug/* +Release/* CMakeUserPresets.json +**/make_patches.sh +.fortls +spack-* diff --git a/.gitmodules b/.gitmodules index 618e7a7..5b2a41f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -16,3 +16,6 @@ [submodule "WW3/WW3"] path = WW3/WW3 url = https://github.com/ESCOMP/WW3 +[submodule "MOM6/GFDL_generic_tracers"] + path = MOM6/GFDL_generic_tracers + url = https://github.com/ACCESS-NRI/GFDL-generic-tracers.git diff --git a/MOM6/CMakeLists.txt b/MOM6/CMakeLists.txt index 8e03a62..a51bbb0 100644 --- a/MOM6/CMakeLists.txt +++ b/MOM6/CMakeLists.txt @@ -7,6 +7,11 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") set(fortran_compile_flags -r8) endif() +add_compile_definitions( + _USE_GENERIC_TRACER + _USE_MOM6_DIAG +) + ### Targets ## MOM6 library @@ -133,7 +138,6 @@ target_sources(OM3_mom6 PRIVATE MOM6/src/framework/MOM_array_transform.F90 MOM6/src/framework/MOM_checksums.F90 MOM6/src/framework/MOM_coms.F90 - MOM6/src/framework/MOM_coupler_types.F90 MOM6/src/framework/MOM_cpu_clock.F90 MOM6/src/framework/MOM_data_override.F90 MOM6/src/framework/MOM_diag_mediator.F90 @@ -238,7 +242,6 @@ target_sources(OM3_mom6 PRIVATE MOM6/src/tracer/ideal_age_example.F90 MOM6/src/tracer/ISOMIP_tracer.F90 MOM6/src/tracer/MOM_CFC_cap.F90 - MOM6/src/tracer/MOM_generic_tracer.F90 MOM6/src/tracer/MOM_hor_bnd_diffusion.F90 MOM6/src/tracer/MOM_neutral_diffusion.F90 MOM6/src/tracer/MOM_OCMIP2_CFC.F90 @@ -298,9 +301,46 @@ target_sources(OM3_mom6 PRIVATE MOM6/config_src/external/drifters/MOM_particles.F90 MOM6/config_src/external/drifters/MOM_particles_types.F90 - MOM6/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 - MOM6/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 - MOM6/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 + + GFDL_generic_tracers/generic_tracers/FMS_coupler_util.F90 + GFDL_generic_tracers/generic_tracers/generic_tracer.F90 + GFDL_generic_tracers/generic_tracers/generic_tracer_utils.F90 + GFDL_generic_tracers/generic_tracers/FMS_ocmip2_co2calc.F90 + GFDL_generic_tracers/generic_tracers/generic_abiotic.F90 + GFDL_generic_tracers/generic_tracers/generic_age.F90 + GFDL_generic_tracers/generic_tracers/generic_argon.F90 + GFDL_generic_tracers/generic_tracers/generic_BLING.F90 + GFDL_generic_tracers/generic_tracers/generic_blres.F90 + GFDL_generic_tracers/generic_tracers/generic_CFC.F90 + GFDL_generic_tracers/generic_tracers/generic_COBALT.F90 + GFDL_generic_tracers/generic_tracers/generic_ERGOM.F90 + GFDL_generic_tracers/generic_tracers/generic_miniBLING.F90 + GFDL_generic_tracers/generic_tracers/generic_SF6.F90 + GFDL_generic_tracers/generic_tracers/generic_TOPAZ.F90 + + GFDL_generic_tracers/mocsy/src/mocsy_buffesm.F90 + GFDL_generic_tracers/mocsy/src/mocsy_constants.F90 + GFDL_generic_tracers/mocsy/src/mocsy_depth2press.F90 + GFDL_generic_tracers/mocsy/src/mocsy_derivauto.F90 + GFDL_generic_tracers/mocsy/src/mocsy_derivnum.F90 + GFDL_generic_tracers/mocsy/src/mocsy_DNAD.F90 + GFDL_generic_tracers/mocsy/src/mocsy_errors.F90 + GFDL_generic_tracers/mocsy/src/mocsy_f2pCO2.F90 + GFDL_generic_tracers/mocsy/src/mocsy_gasx.F90 + GFDL_generic_tracers/mocsy/src/mocsy_p2fCO2.F90 + GFDL_generic_tracers/mocsy/src/mocsy_p80.F90 + GFDL_generic_tracers/mocsy/src/mocsy_phsolvers.F90 + GFDL_generic_tracers/mocsy/src/mocsy_rho.F90 + GFDL_generic_tracers/mocsy/src/mocsy_rhoinsitu.F90 + GFDL_generic_tracers/mocsy/src/mocsy_singledouble.F90 + GFDL_generic_tracers/mocsy/src/mocsy_sw_adtg.F90 + GFDL_generic_tracers/mocsy/src/mocsy_sw_ptmp.F90 + GFDL_generic_tracers/mocsy/src/mocsy_sw_temp.F90 + GFDL_generic_tracers/mocsy/src/mocsy_tis.F90 + GFDL_generic_tracers/mocsy/src/mocsy_tpot.F90 + GFDL_generic_tracers/mocsy/src/mocsy_vars.F90 + GFDL_generic_tracers/mocsy/src/mocsy_varsolver.F90 + MOM6/config_src/external/ODA_hooks/kdtree.f90 MOM6/config_src/external/ODA_hooks/ocean_da_core.F90 MOM6/config_src/external/ODA_hooks/ocean_da_types.F90 @@ -311,7 +351,6 @@ target_sources(OM3_mom6 PRIVATE MOM6/config_src/infra/FMS2/MOM_coms_infra.F90 MOM6/config_src/infra/FMS2/MOM_constants.F90 - MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90 MOM6/config_src/infra/FMS2/MOM_cpu_clock_infra.F90 MOM6/config_src/infra/FMS2/MOM_data_override_infra.F90 MOM6/config_src/infra/FMS2/MOM_diag_manager_infra.F90 @@ -322,14 +361,19 @@ target_sources(OM3_mom6 PRIVATE MOM6/config_src/infra/FMS2/MOM_io_infra.F90 MOM6/config_src/infra/FMS2/MOM_time_manager.F90 - MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90 MOM6/config_src/drivers/nuopc_cap/mom_cap_time.F90 MOM6/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 MOM6/config_src/drivers/nuopc_cap/ocn_comp_NUOPC.F90 MOM6/config_src/drivers/nuopc_cap/time_utils.F90 + + extra_sources/mom_cap_gtracer_flux.F90 ) add_patched_source(OM3_mom6 MOM6/config_src/drivers/nuopc_cap/mom_cap.F90) +add_patched_source(OM3_mom6 MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90) add_patched_source(OM3_mom6 MOM6/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90) +add_patched_source(OM3_mom6 MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90) +add_patched_source(OM3_mom6 MOM6/src/framework/MOM_coupler_types.F90) +add_patched_source(OM3_mom6 MOM6/src/tracer/MOM_generic_tracer.F90) ### Install and Export diff --git a/MOM6/GFDL_generic_tracers b/MOM6/GFDL_generic_tracers new file mode 160000 index 0000000..778750d --- /dev/null +++ b/MOM6/GFDL_generic_tracers @@ -0,0 +1 @@ +Subproject commit 778750dfa69b94c4237b9801a2e7537fea9552ef diff --git a/MOM6/extra_sources/mom_cap_gtracer_flux.F90 b/MOM6/extra_sources/mom_cap_gtracer_flux.F90 new file mode 100644 index 0000000..4dcbd7e --- /dev/null +++ b/MOM6/extra_sources/mom_cap_gtracer_flux.F90 @@ -0,0 +1,674 @@ +!> 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, 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/77618869f48507c8629f28457cb701e25e1ea4fc/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/main/full/coupler_main.F90#L1896 +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/main/full/coupler_main.F90#L2086 +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/77618869f48507c8629f28457cb701e25e1ea4fc/full/coupler_main.F90#L2026 +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/7f585284f1487c0629f8075be350385e6e75dd2f/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/6442d387153064644325c96a5e9e2935139d5e3c/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/MOM6/patches/MOM_coupler_types.F90.patch b/MOM6/patches/MOM_coupler_types.F90.patch new file mode 100644 index 0000000..25a2759 --- /dev/null +++ b/MOM6/patches/MOM_coupler_types.F90.patch @@ -0,0 +1,28 @@ +diff --git a/MOM6/src/framework/MOM_coupler_types.F90 b/MOM6/src/framework/MOM_coupler_types.F90.new +index f87b409..124d786 100644 +--- a/MOM6/src/framework/MOM_coupler_types.F90 ++++ b/MOM6/src/framework/MOM_coupler_types.F90.new +@@ -8,7 +8,10 @@ use MOM_couplertype_infra, only : CT_set_diags, CT_send_data, CT_write_chksums, + 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 + +@@ -22,7 +25,10 @@ public :: atmos_ocn_coupler_flux, coupler_type_data_override + 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/MOM6/patches/MOM_couplertype_infra.F90.patch b/MOM6/patches/MOM_couplertype_infra.F90.patch new file mode 100644 index 0000000..2b7cf4f --- /dev/null +++ b/MOM6/patches/MOM_couplertype_infra.F90.patch @@ -0,0 +1,28 @@ +diff --git a/MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90 b/MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90.new +index 3bcccc1..50fd5b0 100644 +--- a/MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90 ++++ b/MOM6/config_src/infra/FMS2/MOM_couplertype_infra.F90.new +@@ -9,7 +9,10 @@ use coupler_types_mod, only : coupler_type_write_chksums, coupler_type_redistrib + 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 @@ public :: CT_set_diags, CT_send_data, CT_data_override, CT_write_chksums + 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/MOM6/patches/mom_cap.F90.patch b/MOM6/patches/mom_cap.F90.patch index 3db4546..42511ad 100644 --- a/MOM6/patches/mom_cap.F90.patch +++ b/MOM6/patches/mom_cap.F90.patch @@ -1,6 +1,271 @@ ---- MOM6/config_src/drivers/nuopc_cap/mom_cap.F90 2023-08-31 10:49:33.115826000 +1000 -+++ MOM6/config_src/drivers/nuopc_cap/mom_cap.F90.new 2023-08-31 10:50:34.446781000 +1000 -@@ -1672,7 +1672,7 @@ +diff --git a/MOM6/config_src/drivers/nuopc_cap/mom_cap.F90 b/MOM6/config_src/drivers/nuopc_cap/mom_cap.F90.new +index 3574943..48e4532 100644 +--- a/MOM6/config_src/drivers/nuopc_cap/mom_cap.F90 ++++ b/MOM6/config_src/drivers/nuopc_cap/mom_cap.F90.new +@@ -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 +@@ -30,6 +31,7 @@ use MOM_cap_methods, only: med2mod_areacor, state_diagnose + 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 @@ use nuopc_shr_methods, only: get_component_instance + #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 +@@ -96,11 +107,14 @@ implicit none; private + 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 +@@ -408,6 +422,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 +@@ -439,6 +457,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=:), allocatable :: rpointer_filename +@@ -517,6 +536,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, & +@@ -648,13 +669,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) +@@ -721,6 +773,31 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) + Ice_ocean_boundary%hcond = 0.0 + 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 ++ + call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method) + if (use_waves) then + if (wave_method == "EFACTOR") then +@@ -789,6 +866,19 @@ 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) then ++ call fld_list_add(fldsToOcn_num, fldsToOcn, stdname, "will provide") ++ else ++ call ESMF_LogWrite('MOM_cap: no field coupled for generic tracer flux "'//& ++ trim(gas_fluxes%bc(n)%name)//'"', ESMF_LOGMSG_WARNING) ++ endif ++ 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") +@@ -800,6 +890,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) + call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Fioo_q" , "will provide") + call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_bldepth" , "will provide") + ++ ! 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 +@@ -1611,11 +1703,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 +@@ -1627,12 +1722,13 @@ 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(240) :: additional_restart_dir + integer :: num_rest_files + real(8) :: MPI_Wtime, timers + logical :: write_restart +@@ -1673,6 +1769,7 @@ subroutine ModelAdvance(gcomp, rc) + + Time_step_coupled = esmf2fms_time(timeStep) + Time = esmf2fms_time(currTime) ++ Time_import = Time + + !--------------- + ! Apply ocean lag for startup runs: +@@ -1748,8 +1845,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 +@@ -1811,7 +1939,7 @@ subroutine ModelAdvance(gcomp, rc) ! determine restart filename call ESMF_ClockGetNextTime(clock, MyTime, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -9,3 +274,45 @@ if (ChkErr(rc,__LINE__,u_FILE_u)) return if (cesm_coupled) then +@@ -1869,6 +1997,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 +@@ -2149,6 +2285,7 @@ subroutine ocean_model_finalize(gcomp, rc) + integer :: alarmCount + character(len=64) :: timestamp + logical :: write_restart ++ character(240) :: additional_restart_dir + character(len=*),parameter :: subname='(MOM_cap:ocean_model_finalize)' + real(8) :: MPI_Wtime, timefs + +@@ -2182,6 +2319,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/MOM6/patches/mom_cap_methods.F90.patch b/MOM6/patches/mom_cap_methods.F90.patch new file mode 100644 index 0000000..63f58a5 --- /dev/null +++ b/MOM6/patches/mom_cap_methods.F90.patch @@ -0,0 +1,96 @@ +diff --git a/MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90.new +index 125bae5..d78b2f9 100644 +--- a/MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90 ++++ b/MOM6/config_src/drivers/nuopc_cap/mom_cap_methods.F90.new +@@ -20,8 +20,15 @@ use MOM_ocean_model_nuopc, only: ocean_public_type, ocean_state_type + 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,16 @@ 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 +@@ -88,7 +100,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 + +@@ -364,6 +379,46 @@ 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)//& ++ ', coupling field '//trim(stdname), ESMF_LOGMSG_INFO) ++ 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 coupled', ESMF_LOGMSG_INFO) ++ endif ++ enddo ++ deallocate(work) ++ endif ++#endif ++ + end subroutine mom_import + + !> Maps outgoing ocean data to ESMF State diff --git a/MOM6/patches/mom_ocean_model_nuopc.F90.patch b/MOM6/patches/mom_ocean_model_nuopc.F90.patch index ba3e091..273a96a 100644 --- a/MOM6/patches/mom_ocean_model_nuopc.F90.patch +++ b/MOM6/patches/mom_ocean_model_nuopc.F90.patch @@ -1,10 +1,10 @@ diff --git a/MOM6/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/MOM6/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90.new -index 9c81a67..38fface 100644 +index 04b60b0..9855c24 100644 --- a/MOM6/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/MOM6/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90.new @@ -19,6 +19,7 @@ use MOM_coms, only : field_chksum use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf - use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging + 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 @@ -28,7 +28,7 @@ index 9c81a67..38fface 100644 !! 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 -@@ -409,14 +408,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i +@@ -410,14 +409,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, OS%restart_CSp) @@ -40,12 +40,11 @@ index 9c81a67..38fface 100644 - 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) ++ 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. -@@ -531,8 +524,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & +@@ -532,8 +524,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. @@ -55,7 +54,7 @@ index 9c81a67..38fface 100644 weight = 1.0 -@@ -794,7 +786,7 @@ end subroutine ocean_model_end +@@ -795,7 +786,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). @@ -64,7 +63,7 @@ index 9c81a67..38fface 100644 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,16 +818,12 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) +@@ -827,16 +818,12 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) end subroutine ocean_model_save_restart !> Initialize the public ocean type @@ -86,7 +85,7 @@ index 9c81a67..38fface 100644 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 -@@ -847,14 +835,9 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, +@@ -848,14 +835,9 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, ! and have no halos. integer :: isc, iec, jsc, jec @@ -104,7 +103,7 @@ index 9c81a67..38fface 100644 allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), & Ocean_sfc%s_surf (isc:iec,jsc:jec), & -@@ -910,8 +893,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ +@@ -911,8 +893,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)