Merge branch 'user/nnz/accumulate_gtracer_fluxes_2015.06.03' into dev…

…/master

config_src/coupled_driver/ocean_model_MOM.F90

@@ -72,7 +72,7 @@ module ocean_model_mod
 #include <MOM_memory.h>
 
 #ifdef _USE_GENERIC_TRACER
-use MOM_generic_tracer, only : MOM_generic_flux_init
+use MOM_generic_tracer, only : MOM_generic_flux_init,MOM_generic_tracer_fluxes_accumulate 
 #endif
 
 implicit none ; private
@@ -351,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)
@@ -374,10 +374,15 @@ 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
 
@@ -397,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

@@ -1339,17 +1339,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
@@ -1475,8 +1476,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
+  
 end subroutine forcing_accumulate
 
 

src/tracer/MOM_generic_tracer.F90

@@ -44,6 +44,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 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
@@ -77,6 +78,7 @@ module MOM_generic_tracer
   public MOM_generic_tracer_stock
   public MOM_generic_flux_init
   public MOM_generic_tracer_min_max
+  public MOM_generic_tracer_fluxes_accumulate
 
   type, public :: MOM_generic_tracer_CS ; private
      character(len = 200) :: IC_file ! The file in which the generic tracer initial values can
@@ -328,6 +330,12 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, h, OBC, CS, sponge_CSp
     g_tracer=>CS%g_tracer_list
 
     do
+      if(INDEX(CS%IC_file, '_NULL_')) then
+         call MOM_error(WARNING,"The name of the IC_file "//trim(CS%IC_file)//&
+                              " indicates no MOM initialization was asked for the generic tracers."//&
+                              "Bypassing the MOM initialization of ALL generic tracers!")
+         exit
+      endif
       call g_tracer_get_alias(g_tracer,g_tracer_name)
       call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field)
       tr_ptr => tr_field(:,:,:,1)
@@ -336,7 +344,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, h, OBC, CS, sponge_CSp
           .not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then
 
         if (len_trim(CS%IC_file) > 0) then
-          !  Read the tracer concentrations from a netcdf file.
+        !  Read the tracer concentrations from a netcdf file.
           if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, &
                   "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file)
           if (CS%Z_IC_file) then
@@ -351,9 +359,10 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, h, OBC, CS, sponge_CSp
                             "initialized generic tracer "//trim(g_tracer_name)//&
                             " using Generic Tracer File on Z: "//CS%IC_file)
           else
-            ! native grid
+            ! native grid 
             call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
-                  "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//".")
+                  "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//&
+                  " for tracer "//trim(g_tracer_name))
             call read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, domain=G%Domain%mpp_domain)
           endif
         else
@@ -494,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)
+    !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
@@ -579,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">
@@ -806,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
 
@@ -846,6 +854,13 @@ subroutine MOM_generic_flux_init
 
   end subroutine MOM_generic_flux_init
 
+  subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight)
+    type(forcing),         intent(in)    :: flux_tmp
+    real,                  intent(in)    :: weight
+
+   call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight)
+
+  end subroutine MOM_generic_tracer_fluxes_accumulate
 
   subroutine MOM_generic_tracer_get(name,member,array, CS)
     character(len=*),         intent(in)  :: name