From 3b760c0acdf15772bc8d55e6eef34d1ebaf02b9d Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Thu, 23 Jan 2025 10:40:45 +1100 Subject: [PATCH 1/3] Replicate FMS cap initialize_ocean_public_type changes in NUOPC cap The NUOPC cap appears to have been written from the MCT cap, which does not include changes needed to successfully register restarts for FMS coupler type. See NCAR/MOM6@73304eb --- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 47 ++++++------------- 1 file changed, 14 insertions(+), 33 deletions(-) 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 9ac40daaa4..b6cc4de0be 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 @@ -51,8 +52,6 @@ module MOM_ocean_model_nuopc 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 +101,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 @@ -407,14 +406,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. @@ -529,8 +521,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 @@ -792,7 +783,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) @@ -823,16 +814,12 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) 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 @@ -844,14 +831,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), & Ocean_sfc%s_surf (isc:iec,jsc:jec), & @@ -907,8 +889,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, & From ae679238d886a1014554f574e6ebed6d53138111 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Thu, 23 Jan 2025 11:36:38 +1100 Subject: [PATCH 2/3] Allow generic tracers in NUOPC cap These changes provide an initial implementation to allow use of coupled generic tracers with NUOPC-coupled MOM6. This is a bit of a hack as coupling in generic tracers effectively assumes the use of FMScoupler. Broadly, the changes to the NUOPC cap in this commit do the following: - Initialisation phase: Initialise the FMS coupler type data structures used by generic tracers for handling surface fluxes. Register with NUOPC the additional import fields required by the generic tracer package being used. Export of additional fields is not currently handled. - Advance phase: Populate the relevant FMS coupler type data structure with fields received from the atmosphere. Calculate the surface fluxes using a modified version of the routine used by FMScoupler. Restart read/write for the FMS coupler types is also handled. Fields in the FMS coupler types can be overridden using the data_table --- config_src/drivers/nuopc_cap/mom_cap.F90 | 153 +++- .../nuopc_cap/mom_cap_gtracer_flux.F90 | 680 ++++++++++++++++++ .../drivers/nuopc_cap/mom_cap_methods.F90 | 60 +- .../infra/FMS1/MOM_couplertype_infra.F90 | 10 +- .../infra/FMS2/MOM_couplertype_infra.F90 | 10 +- src/framework/MOM_coupler_types.F90 | 10 +- 6 files changed, 912 insertions(+), 11 deletions(-) create mode 100644 config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index a9d4c688f8..1bc4823a20 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 @@ -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 @@ -419,6 +433,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 @@ -450,6 +468,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 @@ -528,6 +547,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, & @@ -659,13 +680,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) @@ -732,6 +784,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 @@ -800,6 +877,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") @@ -811,6 +897,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 @@ -1622,11 +1710,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 @@ -1638,12 +1729,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, write_restartfh @@ -1684,6 +1776,7 @@ subroutine ModelAdvance(gcomp, rc) Time_step_coupled = esmf2fms_time(timeStep) Time = esmf2fms_time(currTime) + Time_import = Time !--------------- ! Apply ocean lag for startup runs: @@ -1759,8 +1852,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 @@ -1827,7 +1951,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 @@ -1885,6 +2009,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 @@ -2124,6 +2256,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 @@ -2157,6 +2290,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 125bae5748..3f81f6aa7e 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 @@ -88,7 +101,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 +380,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 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/framework/MOM_coupler_types.F90 b/src/framework/MOM_coupler_types.F90 index b931a2ddd2..cac9309b5e 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 From 2d5251d8045ce4c6756c99a7a1bcf2dda67aaaa4 Mon Sep 17 00:00:00 2001 From: dougiesquire Date: Thu, 23 Jan 2025 13:36:06 +1100 Subject: [PATCH 3/3] Accumulate fluxes%salt_flux_added and pass to new generic_tracer_update_from_coupler routine This allows virtual flux corrections in generic tracer packages --- config_src/external/GFDL_ocean_BGC/generic_tracer.F90 | 9 +++++++++ src/core/MOM_forcing_type.F90 | 5 +++++ src/tracer/MOM_generic_tracer.F90 | 7 ++++++- 3 files changed, 20 insertions(+), 1 deletion(-) 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/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 998713d1f1..a2efcaac7d 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -2280,6 +2280,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/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 6b10d15e2f..e092aa7d6d 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -21,7 +21,7 @@ module MOM_generic_tracer 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 + 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 @@ -527,6 +527,11 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! 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 !