Use new MARBL interface to initialize Chl

Total chlorophyll is set to 0 initially, and is not updated until after the
first time MOM6 calls get_chl_from_model. This causes two problems:
1. We do not want to compute swpen based on 0 chlorophyll in the first time
step
2. Restart tests to fail, because total chlorophyll was being reset back to 0
in the first time step of the restart

There are two related fixes in this commit:
1. All interior tendency / surface flux outputs from MARBL are now in the
restart file, so the model passes ERS tests
2. If total chlorophyll was not initialized from a restart file, it is computed
based on the initial values of the MARBL tracers (which could be from an IC
file or a restart file)

config_src/external/MARBL/marbl_interface.F90

@@ -36,7 +36,8 @@ module marbl_interface
        contains
         procedure, public  :: put_setting                !< dummy put_setting routine
         procedure, public  :: get_setting                !< dummy get_setting routine
-        procedure, public  :: init                       !< dummy routine
+        procedure, public  :: init                       !< dummy init routine
+        procedure, public  :: compute_totChl             !< dummy routine to compute total Chlorophyll
         procedure, public  :: surface_flux_compute       !< dummy surface flux routine
         procedure, public  :: interior_tendency_compute  !< dummy interior tendency routine
         procedure, public  :: add_output_for_GCM         !< dummy add_output_for_GCM routine
@@ -89,6 +90,15 @@ subroutine init(self,                  &
         call MOM_error(FATAL, error_msg)
     end subroutine init
 
+    !> Dummy version of MARBL's compute_totChl() function
+    subroutine compute_totChl(self)
+
+      class(marbl_interface_class), intent(inout) :: self
+
+      call MOM_error(FATAL, error_msg)
+
+    end subroutine compute_totChl
+
     !> Dummy version of MARBL's surface_flux_compute() function
     subroutine surface_flux_compute(self)
 

src/tracer/MARBL_tracers.F90

@@ -92,8 +92,8 @@ module MARBL_tracers
   character(len=200) :: short_name !< name of variable being saved
   character(len=200) :: file_varname !< name of variable in restart file
   character(len=200) :: units !< variable units
-  real, pointer :: field_2d(:,:) !< memory for 2D field
-  real, pointer :: field_3d(:,:,:) !< memory for 3D field
+  real, pointer :: field_2d(:,:) => NULL() !< memory for 2D field
+  real, pointer :: field_3d(:,:,:) => NULL() !< memory for 3D field
 end type saved_state_for_MARBL_type
 
 !> All calls to MARBL are done via the interface class
@@ -198,12 +198,12 @@ module MARBL_tracers
   real :: DIC_salt_ratio !< ratio to convert salt surface flux to DIC surface flux [conc ppt-1]
   real :: ALK_salt_ratio !< ratio to convert salt surface flux to ALK surface flux [conc ppt-1]
 
-  real, allocatable :: STF(:,:,:)    !< surface fluxes returned from MARBL to use in tracer_vertdiff
-                                     !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1]
-  real, allocatable :: SFO(:,:,:)    !< surface flux output returned from MARBL for use in GCM
-                                     !! e.g. CO2 flux to pass to atmosphere (dims: i, j, num_sfo)
-  real, allocatable :: ITO(:,:,:,:)  !< interior tendency output returned from MARBL for use in GCM
-                                     !! e.g. total chlorophyll to use in shortwave penetration (dims: i, j, k, num_ito)
+  real, allocatable :: STF(:,:,:)          !< surface fluxes returned from MARBL to use in tracer_vertdiff
+                                           !! (dims: i, j, tracer) [conc Z T-1 ~> conc m s-1]
+  real, pointer :: SFO(:,:,:) => NULL()    !< surface flux output returned from MARBL for use in GCM
+                                           !! e.g. CO2 flux to pass to atmosphere (dims: i, j, num_sfo)
+  real, pointer :: ITO(:,:,:,:) => NULL()  !< interior tendency output returned from MARBL for use in GCM
+                                           !! e.g. total chlorophyll to use in shortwave penetration (dims: i, j, k, num_ito)
 
   integer :: u10_sqr_ind   !< index of MARBL forcing field array to copy 10-m wind (squared) into
   integer :: sss_ind       !< index of MARBL forcing field array to copy sea surface salinity into
@@ -413,6 +413,8 @@ subroutine configure_MARBL_tracers(GV, US, param_file, CS)
   ! (4) Request fields needed by MOM6
   CS%sfo_cnt = 0
   CS%ito_cnt = 0
+  CS%flux_co2_ind = -1
+  CS%total_Chl_ind = -1
 
   if (CS%base_bio_on) then
     ! CO2 Flux to the atmosphere
@@ -547,7 +549,8 @@ end subroutine configure_MARBL_tracers
 
 !> This subroutine is used to register tracer fields and subroutines
 !! to be used with MOM.
-function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS, MARBL_computes_chl)
+function register_MARBL_tracers(G, HI, GV, US, param_file, CS, tr_Reg, restart_CS, MARBL_computes_chl)
+  type(ocean_grid_type),      intent(in) :: G    !< The ocean's grid structure
   type(hor_index_type),       intent(in) :: HI   !< A horizontal index type structure.
   type(verticalGrid_type),    intent(in) :: GV   !< The ocean's vertical grid structure
   type(unit_scale_type),      intent(in) :: US   !< A dimensional unit scaling type
@@ -809,6 +812,24 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS,
   call setup_saved_state(MARBL_instances%interior_tendency_saved_state, HI, GV, restart_CS, &
       CS%tracers_may_reinit, CS%interior_tendency_saved_state)
 
+  ! Set up memory for additional output from MARBL and add to restart files
+  allocate(CS%SFO(SZI_(G), SZJ_(G), CS%sfo_cnt), &
+           CS%ITO(SZI_(G), SZJ_(G), SZK_(GV), CS%ito_cnt), &
+           source=0.0)
+
+  do m=1,CS%sfo_cnt
+    write(var_name, "(2A)") 'MARBL_SFO_', &
+                            trim(MARBL_instances%surface_flux_output%outputs_for_GCM(m)%short_name)
+    call register_restart_field(CS%SFO(:,:,m), var_name, .false., restart_CS)
+  enddo
+
+  do m=1,CS%ito_cnt
+    write(var_name, "(2A)") 'MARBL_ITO_', &
+                            trim(MARBL_instances%interior_tendency_output%outputs_for_GCM(m)%short_name)
+    call register_restart_field(CS%ITO(:,:,:,m), var_name, .false., restart_CS)
+  enddo
+
+
   CS%tr_Reg => tr_Reg
   CS%restart_CSp => restart_CS
 
@@ -857,8 +878,6 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag
   ! Allocate memory for surface tracer fluxes
   allocate(CS%STF(SZI_(G), SZJ_(G), CS%ntr), &
            CS%RIV_FLUXES(SZI_(G), SZJ_(G), CS%ntr), &
-           CS%SFO(SZI_(G), SZJ_(G), CS%sfo_cnt), &
-           CS%ITO(SZI_(G), SZJ_(G), SZK_(G), CS%ito_cnt), &
            source=0.0)
 
   ! Allocate memory for d14c forcing
@@ -925,6 +944,7 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag
       diag%axesTL, & ! T=> tracer grid? L => layer center
       day, "Conversion Factor for Bottom Flux -> Tend", "1/m")
 
+  ! Initialize tracers (if they weren't initialized from restart file)
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=name, caller="initialize_MARBL_tracers")
     if ((.not. restart) .or. &
@@ -944,6 +964,30 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag
     endif
   enddo
 
+  ! Initialize total chlorophyll to get SW Pen correct (if it wasn't initialized from restart file)
+  if ((CS%total_Chl_ind > 0) .and. &
+      ((.not. restart) .or. &
+       (.not. query_initialized(CS%ITO(:,:,:,CS%total_Chl_ind), "MARBL_ITO_total_Chl", CS%restart_CSp)))) then
+    ! Three steps per column
+    do j=G%jsc, G%jec
+      do i=G%isc, G%iec
+        ! (i) Copy initial tracers into MARBL structure
+        do m=1,CS%ntr
+          do k=1,GV%ke
+            MARBL_instances%tracers(m,k) = max(CS%tracer_data(m)%tr(i,j,k), 0.)
+          enddo
+        enddo
+        ! (ii) Compute total Chl for the column
+        call MARBL_instances%compute_totChl()
+        ! (iii) Copy total Chl from MARBL data-structure into CS%ITO
+        do k=1,GV%ke
+          CS%ITO(i,j,k,CS%total_Chl_ind) = &
+            MARBL_instances%interior_tendency_output%outputs_for_GCM(CS%total_Chl_ind)%forcing_field_1d(1,k)
+        enddo
+      enddo
+    enddo
+  endif
+
   ! Register diagnostics for river fluxes
   CS%no3_riv_flux = register_diag_field("ocean_model", "NO3_RIV_FLUX", &
       diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
@@ -2116,7 +2160,14 @@ subroutine MARBL_tracers_end(CS)
     if (allocated(CS%qsw_cat_id)) deallocate(CS%qsw_cat_id)
     if (allocated(CS%STF)) deallocate(CS%STF)
     if (allocated(CS%RIV_FLUXES)) deallocate(CS%RIV_FLUXES)
-    if (allocated(CS%SFO)) deallocate(CS%SFO)
+    if (associated(CS%SFO)) then
+      deallocate(CS%SFO)
+      nullify(CS%SFO)
+    endif
+    if (associated(CS%ITO)) then
+      deallocate(CS%ITO)
+      nullify(CS%ITO)
+    endif
     if (allocated(CS%tracer_restoring_ind)) deallocate(CS%tracer_restoring_ind)
     if (allocated(CS%tracer_I_tau_ind)) deallocate(CS%tracer_I_tau_ind)
     if (allocated(CS%fesedflux_in)) deallocate(CS%fesedflux_in)

src/tracer/MOM_tracer_flow_control.F90

@@ -255,7 +255,7 @@ subroutine call_tracer_register(G, GV, US, param_file, CS, tr_Reg, restart_CS)
     register_ideal_age_tracer(G%HI, GV, param_file, CS%ideal_age_tracer_CSp, &
                               tr_Reg, restart_CS)
   if (CS%use_MARBL_tracers) CS%use_MARBL_tracers = &
-    register_MARBL_tracers(G%HI, GV, US, param_file,  CS%MARBL_tracers_CSp, &
+    register_MARBL_tracers(G, G%HI, GV, US, param_file,  CS%MARBL_tracers_CSp, &
                         tr_Reg, restart_CS, CS%get_chl_from_MARBL)
   if (CS%use_regional_dyes) CS%use_regional_dyes = &
     register_dye_tracer(G%HI, GV, US, param_file, CS%dye_tracer_CSp, &