Implement flux accumulation for generic tracers

- generic_tracer_coupler_get has to be called from step_MOM() rather than MOM_generic_tracer_column_physics()
  so that the generic tracer fluxes can be accumulated on steps that do not "modify" tracers, particularly
  when tracer timestep spans dt_cpld ( dt_therm > dt_cpld)

- the fluxes need to be accumulated with the same weight that MOM tracers are being accumulated, hence the need
  to expose the "weight" from MOM_forcing_type.F90

- the g_tracer_send_diag calls which writes diagnostics for the generic tracers needs to be relocated to
  where it is called every tracer timestep (and not every dt_cpld) to avoid having wrong instant or mean values
  in their diagnostics, particularly when tracer timestep spans dt_cpld ( dt_therm > dt_cpld)

config_src/coupled_driver/ocean_model_MOM.F90

@@ -72,8 +72,7 @@ module ocean_model_mod
 #include <MOM_memory.h>
 
 #ifdef _USE_GENERIC_TRACER
-use MOM_generic_tracer, only : MOM_generic_flux_init
-use generic_tracer, only: generic_tracer_coupler_get
+use MOM_generic_tracer, only : MOM_generic_flux_init,MOM_generic_tracer_fluxes_accumulate 
 #endif
 
 implicit none ; private
@@ -352,10 +351,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
                                  ! start of a call to step_MOM.
   integer :: index_bnds(4)       ! The computational domain index bounds in the
                                  ! ice-ocean boundary type.
-
+  real :: weight            ! Flux accumulation weight
   real :: time_step         ! The time step of a call to step_MOM in seconds.
   integer :: secs, days
-
+  
   call callTree_enter("update_ocean_model(), ocean_model_MOM.F90")
   call get_time(Ocean_coupling_time_step, secs, days)
   time_step = 86400.0*real(days) + real(secs)
@@ -375,22 +374,22 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), &
                               index_bnds(3), index_bnds(4))
 
+  weight = 1.0
 
   if (OS%fluxes%fluxes_used) then
     call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) ! Needed to allow diagnostics in convert_IOB
     call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, &
                                OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity)
+#ifdef _USE_GENERIC_TRACER
+    call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes
+#endif
 
   ! Add ice shelf fluxes
 
     if (OS%use_ice_shelf) then
       call shelf_calc_flux(OS%State, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp)
     endif
 
-#ifdef _USE_GENERIC_TRACER
-    call generic_tracer_coupler_get(OS%fluxes%tr_fluxes)
-#endif
-
     ! Indicate that there are new unused fluxes.
     OS%fluxes%fluxes_used = .false.
     OS%fluxes%dt_buoy_accum = time_step
@@ -403,7 +402,10 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
       call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp)
     endif
 
-    call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid)
+    call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight)
+#ifdef _USE_GENERIC_TRACER
+    call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average
+#endif
   endif
 
   call disable_averaging(OS%MOM_CSp%diag)

src/core/MOM_forcing_type.F90

@@ -32,9 +32,6 @@ module MOM_forcing_type
 use MOM_variables,     only : surface, thermo_var_ptrs
 
 use coupler_types_mod, only : coupler_2d_bc_type
-#ifdef _USE_GENERIC_TRACER
-  use generic_tracer, only: generic_tracer_coupler_accumulate
-#endif
 
 implicit none ; private
 
@@ -1341,17 +1338,18 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles)
 end subroutine register_forcing_type_diags
 
 
-subroutine forcing_accumulate(flux_tmp, fluxes, dt, G)
+subroutine forcing_accumulate(flux_tmp, fluxes, dt, G, wt2)
   type(forcing),         intent(in)    :: flux_tmp
   type(forcing),         intent(inout) :: fluxes
   real,                  intent(in)    :: dt
   type(ocean_grid_type), intent(inout) :: G
+  real,                  intent(out)   :: wt2
   !   This subroutine copies mechancal forcing from flux_tmp to fluxes and
   ! stores the time-weighted averages of the various buoyancy fluxes in fluxes,
   ! and increments the amount of time over which the buoyancy forcing should be
   ! applied.
 
-  real :: wt1, wt2
+  real :: wt1
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
   is   = G%isc   ; ie   = G%iec    ; js   = G%jsc   ; je   = G%jec
@@ -1477,15 +1475,8 @@ subroutine forcing_accumulate(flux_tmp, fluxes, dt, G)
   endif
 
   !### This needs to be replaced with an appropriate copy and average.
-   fluxes%tr_fluxes => flux_tmp%tr_fluxes
+  fluxes%tr_fluxes => flux_tmp%tr_fluxes
 
-#ifdef _USE_GENERIC_TRACER
-   !Cannot call MOM_generic_tracer_fluxes_accumulate(wt1, flux_tmp%tr_fluxes) since that model "use"s this one
-   !So, call the geneirc_tracer interface directly.
-   !if (CS%use_MOM_generic_tracer)  !CS is not defined here. 
-   call generic_tracer_coupler_accumulate(wt1, flux_tmp%tr_fluxes)
-#endif
-
 end subroutine forcing_accumulate
 
 

src/tracer/MOM_generic_tracer.F90

@@ -503,22 +503,14 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
     !Get the tracer list
     if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL,&
          trim(sub_name)//": No tracer in the list.")
-    !
-    !Update the fields of the generic tracers from Tr(n)
-    !
-    !nnz: This step is not necessary in MOM code because the tracers
-    !are "registered" for restart and advection. I.e., the ocean model
-    !has the pointers to %field arrays and use them directly.
-    !I think this is a much better design than T_prog(n) of MOM which
-    !get updated by ocean model and have to be copied back and for
-    !to the tracers fields.
-    !The problem with MOM is T_prog%field is not a pointer but is an allocatable array.
-
 
     !
     !Extract the tracer surface fields from coupler and update tracer fields from sources
     !
-    !call generic_tracer_coupler_get(fluxes%tr_fluxes) !This is moved out to ocean_model_MOM.F90
+    !call generic_tracer_coupler_get(fluxes%tr_fluxes) 
+    !Niki: This is moved out to ocean_model_MOM.F90 because if dt_therm>dt_cpld we need to average
+    !      the fluxes without coming into this subroutine.
+    !      MOM5 has to modified to conform.
 
     !
     !Add contribution of river to surface flux
@@ -588,6 +580,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
 
     call generic_tracer_update_from_bottom(dt, 1, get_diag_time_end(CS%diag)) !Second arg is tau which is always 1 for MOM
 
+    !Output diagnostics via diag_manager for all generic tracers and their fluxes
+    call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
+
   end subroutine MOM_generic_tracer_column_physics
 
   ! <SUBROUTINE NAME="MOM_generic_tracer_stock">
@@ -815,9 +810,13 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS)
          tau=1)
 
     !Output diagnostics via diag_manager for all tracers in this module
-    if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&
-         "No tracer in the list.")
-    call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
+!    if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&
+!         "No tracer in the list.")
+!    call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
+    !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld
+    !      hence if dt_therm > dt_cpld we get output (and contribution to the mean) at times that tracers 
+    !      had not been updated.
+    !      Moving this to the end of column physics subrotuine fixes this issue.
 
   end subroutine MOM_generic_tracer_surface_state
 
@@ -855,11 +854,11 @@ subroutine MOM_generic_flux_init
 
   end subroutine MOM_generic_flux_init
 
-  subroutine MOM_generic_tracer_fluxes_accumulate(weight, flux_tmp)
-    real,                  intent(in)    :: weight
+  subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight)
     type(forcing),         intent(in)    :: flux_tmp
+    real,                  intent(in)    :: weight
 
-   call generic_tracer_coupler_accumulate(weight, flux_tmp%tr_fluxes)
+   call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight)
 
   end subroutine MOM_generic_tracer_fluxes_accumulate