Add MARBL to MOM6

.gitignore

@@ -4,6 +4,15 @@
 html
 
 
+# Build output
+*.o
+*.mod
+MOM6
+build/
+deps/
+pkg/MARBL
+
+
 # Autoconf output
 aclocal.m4
 autom4te.cache/

config_src/drivers/mct_cap/mom_surface_forcing_mct.F90

@@ -38,6 +38,8 @@ module MOM_surface_forcing_mct
 use time_interp_external_mod, only : init_external_field, time_interp_external
 use time_interp_external_mod, only : time_interp_external_init
 use MOM_io,               only: stdout
+use marbl_forcing_type_mod, only : marbl_forcing_CS, marbl_forcing_init, marbl_forcing_type_init
+use marbl_forcing_type_mod, only : marbl_ice_ocean_boundary_type, convert_marbl_IOB_to_forcings
 
 implicit none ; private
 
@@ -142,6 +144,8 @@ module MOM_surface_forcing_mct
   type(forcing_diags), public :: handles !< diagnostics handles
   type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< restart pointer
   type(user_revise_forcing_CS), pointer :: urf_CS => NULL() !< user revise pointer
+
+  type(marbl_forcing_CS), pointer :: marbl_forcing_CSp => NULL() !< parameters for getting MARBL forcing
 end type surface_forcing_CS
 
 !> Structure corresponding to forcing, but with the elements, units, and conventions
@@ -185,6 +189,9 @@ module MOM_surface_forcing_mct
                                                               !! flux-exchange code, based on what the sea-ice
                                                               !! model is providing.  Otherwise, the value from
                                                               !! the surface_forcing_CS is used.
+
+  type(marbl_ice_ocean_boundary_type), pointer :: MARBL_IOB => NULL() !< Structure containing IOB fields
+                                                                      !! (only needed by MARBL)
 end type ice_ocean_boundary_type
 
 integer :: id_clock_forcing
@@ -313,6 +320,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
     if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
 
+    call marbl_forcing_type_init(isd,ied,jsd,jed,fluxes%MARBL_forcing)
+
   endif   ! endif for allocation and initialization
 
 
@@ -524,6 +533,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
   enddo; enddo
 
+  ! Copy MARBL-specific IOB fields into fluxes%MARBL_forcing
+  call convert_marbl_IOB_to_forcings(IOB%MARBL_IOB, Time, G, US, i0, j0, fluxes%MARBL_forcing, CS%marbl_forcing_CSp)
+
   ! applied surface pressure from atmosphere and cryosphere
   if (associated(IOB%p)) then
      if (CS%max_p_surf >= 0.0) then
@@ -1295,6 +1307,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
     call data_override_init(Ocean_domain_in=G%Domain%mpp_domain)
   endif
 
+  ! Set up MARBL forcing control structure
+  call marbl_forcing_init(G, param_file, CS%diag, Time, CS%inputdir, CS%marbl_forcing_CSp)
+
   if (present(restore_salt)) then ; if (restore_salt) then
     salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file)
     CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain)

config_src/drivers/mct_cap/ocn_cap_methods.F90

@@ -21,7 +21,7 @@ module ocn_cap_methods
 contains
 !=======================================================================
 
-!> Maps incomping ocean data to MOM6 data structures
+!> Maps incoming ocean data to MOM6 data structures
 subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit, Eclock, c1, c2, c3, c4)
   real(kind=8)                  , intent(in)    :: x2o(:,:)           !< incoming data
   type(cpl_indices_type)        , intent(in)    :: ind                !< Structure with MCT attribute vects and indices
@@ -105,6 +105,35 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
         ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(i,j)
         ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(i,j)
       endif
+
+      ! MARBL-specific forcing fields
+      ! Dust and iron fluxes (Do I need the GRID%mask2dT terms? I have them in convert_IOB_to_fluxes())
+      ! TODO: these are in units kg/m^2/s, want g/cm^2/s => multiply by 0.1 somewhere
+      !       currently the conversion is done in MARBL_tracers:MARBL_tracers_column_physics when copying
+      !       into surface_flux_forcings but maybe we should convert somewhere else because we want to
+      ! register these IOB fields as diagnostics (in g/cm^2/s for now)
+      ice_ocean_boundary%MARBL_IOB%atm_fine_dust_flux(i,j) = (x2o(ind%x2o_Faxa_dstwet1,k) + &
+                                                              x2o(ind%x2o_Faxa_dstdry1,k)) * &
+                                                             GRID%mask2dT(i,j)
+      ice_ocean_boundary%MARBL_IOB%atm_coarse_dust_flux(i,j) = (x2o(ind%x2o_Faxa_dstwet2,k) + &
+                                                                x2o(ind%x2o_Faxa_dstwet3,k) + &
+                                                                x2o(ind%x2o_Faxa_dstwet4,k) + &
+                                                                x2o(ind%x2o_Faxa_dstdry2,k) + &
+                                                                x2o(ind%x2o_Faxa_dstdry3,k) + &
+                                                                x2o(ind%x2o_Faxa_dstdry4,k)) * &
+                                                               GRID%mask2dT(i,j)
+      ice_ocean_boundary%MARBL_IOB%seaice_dust_flux(i,j) = x2o(ind%x2o_Fioi_flxdst,k) * GRID%mask2dT(i,j)
+      ice_ocean_boundary%MARBL_IOB%atm_bc_flux(i,j) = (x2o(ind%x2o_Faxa_bcphidry,k) + x2o(ind%x2o_Faxa_bcphodry,k) + &
+                                                       x2o(ind%x2o_Faxa_bcphiwet,k)) * GRID%mask2dT(i,j)
+      ice_ocean_boundary%MARBL_IOB%seaice_bc_flux(i,j) = (x2o(ind%x2o_Fioi_bcpho,k) + x2o(ind%x2o_Fioi_bcphi,k)) * &
+                                                         GRID%mask2dT(i,j)
+
+      ! ice fraction
+      ice_ocean_boundary%MARBL_IOB%ice_fraction(i,j) = x2o(ind%x2o_Si_ifrac,k) * GRID%mask2dT(i,j)
+
+      ! 10m wind
+      ice_ocean_boundary%MARBL_IOB%u10_sqr(i,j) = x2o(ind%x2o_So_duu10n,k) * GRID%mask2dT(i,j)
+
     enddo
   enddo
 

config_src/drivers/mct_cap/ocn_comp_mct.F90

@@ -44,6 +44,8 @@ module ocn_comp_mct
 use mpp_domains_mod,      only: mpp_get_compute_domain
 use MOM_io,               only: stdout
 
+use marbl_forcing_type_mod, only: marbl_iob_allocate
+
 ! Previously inlined - now in separate modules
 use MOM_ocean_model_mct,     only: ocean_public_type, ocean_state_type
 use MOM_ocean_model_mct,     only: ocean_model_init , update_ocean_model, ocean_model_end
@@ -812,8 +814,6 @@ end subroutine ocean_model_init_sfc
 !! mi, mass of ice (kg/m2)
 !!
 !! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in fluxes):
-!! x2o_Si_ifrac, fractional ice wrt ocean
-!! x2o_So_duu10n, 10m wind speed squared (m^2/s^2)
 !! x2o_Sa_co2prog, bottom atm level prognostic CO2
 !! x2o_Sa_co2diag, bottom atm level diagnostic CO2
 !!
@@ -907,6 +907,8 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
   IOB%mi               = 0.0
   IOB%p                = 0.0
 
+  call marbl_iob_allocate(isc, iec, jsc, jec, IOB%MARBL_IOB)
+
 end subroutine IOB_allocate
 
 end module ocn_comp_mct

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -44,6 +44,7 @@ module MOM_cap_mod
 use shr_mpi_mod,              only : shr_mpi_min, shr_mpi_max
 #endif
 use time_utils_mod,           only: esmf2fms_time
+use marbl_forcing_type_mod,   only: marbl_iob_allocate
 
 use, intrinsic :: iso_fortran_env, only: output_unit
 
@@ -687,6 +688,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     Ice_ocean_boundary%vstkb           = 0.0
   endif
 
+  call marbl_iob_allocate(isc, iec, jsc, jec, Ice_ocean_boundary%MARBL_IOB)
+
   ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state
   call ESMF_GridCompSetInternalState(gcomp, ocean_internalstate, rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -727,6 +730,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"                  , "will provide") !-> ice runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn"        , "will provide")
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac"                   , "will provide") !-> ice runoff
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"                  , "will provide") !-> ice runoff
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_flxdst"                , "will provide") !-> ice runoff
   !These are not currently used and changing requires a nuopc dictionary change
   !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx"        , "will provide")
   !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx"       , "will provide")

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -305,6 +305,43 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
     deallocate(stkx1,stkx2,stkx3,stky1,stky2,stky3)
   endif
 
+  ! MARBL fields
+
+  !----
+  ! ice fraction
+  !----
+  ice_ocean_boundary%MARBL_IOB%ice_fraction(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'Si_ifrac',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%MARBL_IOB%ice_fraction, rc=rc)
+  if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+       line=__LINE__, &
+       file=__FILE__)) &
+       return  ! bail out
+
+  !----
+  ! 10m wind squared
+  !----
+  ice_ocean_boundary%MARBL_IOB%u10_sqr(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'So_duu10n',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%MARBL_IOB%u10_sqr, rc=rc)
+  if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+       line=__LINE__, &
+       file=__FILE__)) &
+       return  ! bail out
+
+  !----
+  ! seaice_dust_flux is a single field from the coupler
+  ! atm_fine_dust_flux, atm_coarse_dust_flux, atm_bc_flux, and seaice_bc_flux
+  ! are all sums of multiple fields and will be treated slightly differently
+  !----
+  ice_ocean_boundary%MARBL_IOB%seaice_dust_flux(:,:) = 0._ESMF_KIND_R8
+  call state_getimport(importState, 'Fioi_flxdst',  &
+       isc, iec, jsc, jec, ice_ocean_boundary%MARBL_IOB%seaice_dust_flux, rc=rc)
+  if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+       line=__LINE__, &
+       file=__FILE__)) &
+       return  ! bail out
+
 end subroutine mom_import
 
 !> Maps outgoing ocean data to ESMF State

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -39,6 +39,8 @@ module MOM_surface_forcing_nuopc
 use mpp_mod,              only : mpp_chksum
 use time_interp_external_mod, only : init_external_field, time_interp_external
 use time_interp_external_mod, only : time_interp_external_init
+use marbl_forcing_type_mod,   only : marbl_forcing_CS, marbl_forcing_init, marbl_forcing_type_init
+use marbl_forcing_type_mod,   only : marbl_ice_ocean_boundary_type, convert_marbl_IOB_to_forcings
 
 implicit none ; private
 
@@ -150,6 +152,8 @@ module MOM_surface_forcing_nuopc
 
   type(MOM_restart_CS), pointer :: restart_CSp => NULL()
   type(user_revise_forcing_CS), pointer :: urf_CS => NULL()
+
+  type(marbl_forcing_CS), pointer :: marbl_forcing_CSp => NULL() !< parameters for getting MARBL forcing
 end type surface_forcing_CS
 
 !> Structure corresponding to forcing, but with the elements, units, and conventions
@@ -201,6 +205,9 @@ module MOM_surface_forcing_nuopc
                                                               !! flux-exchange code, based on what the sea-ice
                                                               !! model is providing.  Otherwise, the value from
                                                               !! the surface_forcing_CS is used.
+
+  type(marbl_ice_ocean_boundary_type), pointer :: MARBL_IOB => NULL() !< Structure containing IOB fields
+                                                                      !! (only needed by MARBL)
 end type ice_ocean_boundary_type
 
 integer :: id_clock_forcing
@@ -327,6 +334,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
     if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed)
 
+    call marbl_forcing_type_init(isd,ied,jsd,jed,fluxes%MARBL_forcing)
+
   endif   ! endif for allocation and initialization
 
   if (((associated(IOB%ustar_berg) .and. (.not.associated(fluxes%ustar_berg))) &
@@ -534,6 +543,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
   enddo ; enddo
 
+  ! Copy MARBL-specific IOB fields into fluxes%MARBL_forcing
+  call convert_marbl_IOB_to_forcings(IOB%MARBL_IOB, Time, G, US, i0, j0, fluxes%MARBL_forcing, CS%marbl_forcing_CSp)
+
   ! applied surface pressure from atmosphere and cryosphere
   if (associated(IOB%p)) then
      if (CS%max_p_surf >= 0.0) then
@@ -1335,6 +1347,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
     call data_override_init(Ocean_domain_in=G%Domain%mpp_domain)
   endif
 
+  ! Set up MARBL forcing control structure
+  call marbl_forcing_init(G, param_file, CS%diag, Time, CS%inputdir, CS%marbl_forcing_CSp)
+
   if (present(restore_salt)) then ; if (restore_salt) then
     salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file)
     CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain)

config_src/external/MARBL/README.md

@@ -0,0 +1,6 @@
+GFDL_ocean_BGC
+==============
+
+These APIs reflect those for the MARBL library available at https://github.com/marbl-ecosys/MARBL
+
+The modules in this directory do not do any computations. They simply reflect the APIs of the above package.

config_src/external/MARBL/marbl_interface.F90

@@ -0,0 +1,102 @@
+!> A non-functioning template of the MARBL interface
+module marbl_interface
+
+    use MOM_error_handler,            only : MOM_error, FATAL
+    use marbl_logging,                only : marbl_log_type
+    use marbl_interface_public_types, only : marbl_forcing_fields_type
+    use marbl_interface_public_types, only : marbl_tracer_metadata_type
+    use marbl_interface_public_types, only : marbl_saved_state_type
+    use marbl_interface_public_types, only : marbl_diagnostics_type
+    use marbl_interface_public_types, only : marbl_domain_type
+    implicit none
+    private ! Only want marbl_interface_class to be public, not supporting functions
+
+    !> A non-functioning template of the MARBL_interface class
+    !!
+    !> All variables are dummy representations of actual members of the real marbl_interface_class
+    !! that are used in the MARBL tracer routines.
+    type, public :: marbl_interface_class
+        type(marbl_log_type) :: StatusLog  !< dummy log
+        type(marbl_forcing_fields_type), allocatable :: surface_flux_forcings(:)  !< dummy forcing array
+        type(marbl_forcing_fields_type), allocatable :: interior_tendency_forcings(:)  !< dummy forcing array
+        type(marbl_tracer_metadata_type), allocatable :: tracer_metadata(:)  !< dummy metadata array
+        type(marbl_domain_type) :: domain  !< dummy domain
+        type(marbl_saved_state_type) :: surface_flux_saved_state  !< dummy saved state
+        type(marbl_saved_state_type) :: interior_tendency_saved_state  !< dummy saved state
+        type(marbl_diagnostics_type) :: interior_tendency_diags  !< dummy diagnostics
+        type(marbl_diagnostics_type) :: surface_flux_diags  !< dummy diagnostics
+        real, allocatable :: tracers(:,:)  !< dummy tracer array
+        real, allocatable :: tracers_at_surface(:,:)  !< dummy tracer surface array
+        real, allocatable :: surface_fluxes(:,:)  !< dummy fluxes
+        real, allocatable :: interior_tendencies(:,:)  !< dummy tendencies
+       contains
+        procedure, public  :: init  !< dummy routine
+        procedure, public  :: surface_flux_compute  !< dummy surface flux routine
+        procedure, public  :: interior_tendency_compute  !< dummy interior tendency routine
+        procedure, public  :: shutdown  !< dummy shutdown routine
+        procedure, public :: put_setting  !< dummy put_setting routine
+    end type marbl_interface_class
+
+    !> Error message that appears if the dummy interface is called
+    character(len=*), parameter :: error_msg = "MOM6 built the MARBL stubs rather than the full library"
+
+contains
+
+    !> Dummy version of MARBL's put_setting() function
+    subroutine put_setting(self, str_in)
+        class(marbl_interface_class), intent(in) :: self
+        character(len=*),             intent(in) :: str_in
+
+        call MOM_error(FATAL, error_msg)
+    end subroutine put_setting
+
+    !> Dummy version of MARBL's init() function
+    subroutine init(self,                  &
+        gcm_num_levels,                    &
+        gcm_num_PAR_subcols,               &
+        gcm_num_elements_surface_flux,     &
+        gcm_delta_z,                       &
+        gcm_zw,                            &
+        gcm_zt,                            &
+        lgcm_has_global_ops)
+
+        class(marbl_interface_class), intent(inout) :: self
+        integer,                      intent(in)    :: gcm_num_levels
+        integer,                      intent(in)    :: gcm_num_PAR_subcols
+        integer,                      intent(in)    :: gcm_num_elements_surface_flux
+        real,                         intent(in)    :: gcm_delta_z(gcm_num_levels)
+        real,                         intent(in)    :: gcm_zw(gcm_num_levels)
+        real,                         intent(in)    :: gcm_zt(gcm_num_levels)
+        logical,                      intent(in)    :: lgcm_has_global_ops
+
+        call MOM_error(FATAL, error_msg)
+    end subroutine init
+
+    !> Dummy version of MARBL's surface_flux_compute() function
+    subroutine surface_flux_compute(self)
+
+        class(marbl_interface_class), intent(inout) :: self
+
+        call MOM_error(FATAL, error_msg)
+
+    end subroutine surface_flux_compute
+
+    !> Dummy version of MARBL's interior_tendency_compute() function
+    subroutine interior_tendency_compute(self)
+
+        class(marbl_interface_class), intent(inout) :: self
+
+        call MOM_error(FATAL, error_msg)
+
+    end subroutine interior_tendency_compute
+
+    !> Dummy version of MARBL's shutdown() function
+    subroutine shutdown(self)
+
+        class(marbl_interface_class), intent(inout) :: self
+
+        call MOM_error(FATAL, error_msg)
+
+    end subroutine shutdown
+
+end module marbl_interface