Add generic tracer budget diagnostics (mom-ocean#12)

* Apply changes for BGC budgets

# Conflicts:
#	src/tracer/MOM_generic_tracer.F90
#	src/tracer/MOM_tracer_registry.F90

* Updated budget code from Fan

* Enable built in diagnostics

* Cleanup

* Cleanup

* Add boundary_forcing_tend to stub g_tracer_type

* Clean trailing whitespace

* Remove redundant zeroing of boundary_forcing_tend

* Add an attribution

* Clarify which diagnostics are for concentration

* Fix units for _diffusion_xy

config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90

@@ -20,13 +20,15 @@ module g_tracer_utils
     real, pointer, dimension(:,:,:,:) :: field  => NULL()
     !> Tracer concentration in river runoff
     real, allocatable, dimension(:,:) :: trunoff
+    real, allocatable, dimension(:,:,:)  :: boundary_forcing_tend !< Tendency for budget diagnostics
     logical :: requires_restart = .true. !< Unknown
     character(len=fm_string_len) :: src_file !< Tracer source filename
     character(len=fm_string_len) :: src_var_name !< Tracer source variable name
     character(len=fm_string_len) :: src_var_unit !< Tracer source variable units
     character(len=fm_string_len) :: src_var_gridspec !< Tracer source grid file name
     character(len=fm_string_len) :: obc_src_file_name !< Boundary condition tracer source filename
     character(len=fm_string_len) :: obc_src_field_name !< Boundary condition tracer source fieldname
+    integer :: diag_id_boundary_forcing_tend = -1 !< Budget diagnostic id
     integer :: src_var_record !< Unknown
     logical :: runoff_added_to_stf = .false. !< Has flux in from runoff been added to stf?
     logical :: requires_src_info = .false. !< Unknown

src/tracer/MOM_generic_tracer.F90

@@ -219,7 +219,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
       if (g_tracer_is_prog(g_tracer)) then
         call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, &
                              name=g_tracer_name, longname=longname, units=units, &
-                             registry_diags=.false., &   !### CHANGE TO TRUE?
+                             registry_diags=.true., &   !### CHANGE TO TRUE?
                              restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit)
       else
         call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, &
@@ -673,13 +673,22 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
         if (g_tracer_is_prog(g_tracer)) then
           do k=1,nk ;do j=jsc,jec ; do i=isc,iec
             h_work(i,j,k) = h_old(i,j,k)
+            if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then
+               g_tracer%boundary_forcing_tend(i,j,k) = g_tracer%field(i,j,k,1)
+            endif
           enddo ; enddo ; enddo
           call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, &
                             fluxes, h_work, evap_CFL_limit, minimum_forcing_depth)
+          if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then
+             do k=1,nk ;do j=jsc,jec ; do i=isc,iec
+                g_tracer%boundary_forcing_tend(i,j,k)=G%mask2dT(i,j)*(g_tracer%field(i,j,k,1) &
+                 - g_tracer%boundary_forcing_tend(i,j,k))/dt
+             enddo ; enddo ; enddo
+          endif
         endif
 
-         !traverse the linked list till hit NULL
-         call g_tracer_get_next(g_tracer, g_tracer_next)
+        !traverse the linked list till hit NULL
+        call g_tracer_get_next(g_tracer, g_tracer_next)
         if (.NOT. associated(g_tracer_next)) exit
         g_tracer=>g_tracer_next
       enddo

src/tracer/MOM_tracer_advect.F90

@@ -193,6 +193,9 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
     if (associated(Reg%Tr(m)%ad_x)) Reg%Tr(m)%ad_x(:,:,:) = 0.0
     if (associated(Reg%Tr(m)%ad_y)) Reg%Tr(m)%ad_y(:,:,:) = 0.0
     if (associated(Reg%Tr(m)%advection_xy)) Reg%Tr(m)%advection_xy(:,:,:) = 0.0
+    if (associated(Reg%Tr(m)%advectionc_xy)) Reg%Tr(m)%advectionc_xy(:,:,:) = 0.0
+    if (associated(Reg%Tr(m)%advectionc_x)) Reg%Tr(m)%advectionc_x(:,:,:) = 0.0
+    if (associated(Reg%Tr(m)%advectionc_y)) Reg%Tr(m)%advectionc_y(:,:,:) = 0.0
     if (associated(Reg%Tr(m)%ad2d_x)) Reg%Tr(m)%ad2d_x(:,:) = 0.0
     if (associated(Reg%Tr(m)%ad2d_y)) Reg%Tr(m)%ad2d_y(:,:) = 0.0
   enddo
@@ -365,6 +368,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
                         ! the grid box, both in [H L2 ~> m3 or kg].
   real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the
                         ! current iteration [H L2 ~> m3 or kg].
+  real, dimension(SZIB_(G)) :: tprev !< tracer conc at the end of previous step.
   real, dimension(SZIB_(G)) :: &
     hlst, &             ! Work variable [H L2 ~> m3 or kg].
     Ihnew, &            ! Work variable [H-1 L-2 ~> m-3 or kg-1].
@@ -655,6 +659,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
       do i=is,ie
         if (do_i(i,j)) then
           if (Ihnew(i) > 0.0) then
+            tprev(i)=Tr(m)%t(i,j,k)
             Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - &
                               (flux_x(I,j,m) - flux_x(I-1,j,m))) * Ihnew(i)
           endif
@@ -674,7 +679,16 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
                                           Idt * G%IareaT(i,j)
         endif ; enddo
       endif
-
+      if (associated(Tr(m)%advectionc_xy)) then
+        do i=is,ie ; if (do_i(i,j)) then
+          Tr(m)%advectionc_xy(i,j,k) = Tr(m)%advectionc_xy(i,j,k)+(Tr(m)%t(i,j,k) - tprev(i))*Idt*G%mask2dT(i,j)
+        endif ; enddo
+      endif
+      if (associated(Tr(m)%advectionc_x)) then
+        do i=is,ie ; if (do_i(i,j)) then
+          Tr(m)%advectionc_x(i,j,k) =(Tr(m)%t(i,j,k) - tprev(i))*Idt*G%mask2dT(i,j)
+        endif ; enddo
+      endif
     enddo
 
   endif ; enddo ! End of j-loop.
@@ -736,6 +750,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
     T_tmp               ! The copy of the tracer concentration at constant i,k [conc].
   real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
                                 ! current iteration [H L2 ~> m3 or kg].
+  real, dimension(SZIB_(G)) :: tprev !< tracer conc at the end of previous step.
   real :: hup, hlos             ! hup is the upwind volume, hlos is the
                                 ! part of that volume that might be lost
                                 ! due to advection out the other side of
@@ -1042,10 +1057,10 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
     ! update tracer and save some diagnostics
     do m=1,ntr
       do i=is,ie ; if (do_i(i,j)) then
+        tprev(i)=Tr(m)%t(i,j,k)
         Tr(m)%t(i,j,k) = (Tr(m)%t(i,j,k) * hlst(i) - &
                           (flux_y(i,m,J) - flux_y(i,m,J-1))) * Ihnew(i)
       endif ; enddo
-
       ! diagnose convergence of flux_y and add to convergence of flux_x.
       ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
       if (associated(Tr(m)%advection_xy)) then
@@ -1054,7 +1069,16 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
                                           G%IareaT(i,j)
         endif ; enddo
       endif
-
+      if (associated(Tr(m)%advectionc_xy)) then
+        do i=is,ie ; if (do_i(i,j)) then
+          Tr(m)%advectionc_xy(i,j,k) = Tr(m)%advectionc_xy(i,j,k)+(Tr(m)%t(i,j,k) - tprev(i))*Idt*G%mask2dT(i,j)
+        endif ; enddo
+      endif
+      if (associated(Tr(m)%advectionc_y)) then
+        do i=is,ie ; if (do_i(i,j)) then
+          Tr(m)%advectionc_y(i,j,k) = (Tr(m)%t(i,j,k) - tprev(i))*Idt*G%mask2dT(i,j)
+        endif ; enddo
+      endif
     enddo
   endif ; enddo ! End of j-loop.
 

src/tracer/MOM_tracer_hor_diff.F90

@@ -400,6 +400,16 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
     if (associated(Reg%Tr(m)%df2d_y)) then
       do J=js-1,je ; do i=is,ie ; Reg%Tr(m)%df2d_y(i,J) = 0.0 ; enddo ; enddo
     endif
+    if (associated(Reg%Tr(m)%diffusionc_xy)) then
+      do k=1,nz ; do j=js,je ; do i=is,ie
+        Reg%Tr(m)%diffusionc_xy(i,j,k) = 0.0
+      enddo ; enddo ; enddo
+    endif
+    if (associated(Reg%Tr(m)%diffusion_xy)) then
+      do k=1,nz ; do j=js,je ; do i=is,ie
+        Reg%Tr(m)%diffusion_xy(i,j,k) = 0.0
+      enddo ; enddo ; enddo
+    endif
   enddo
 
   if (CS%use_hor_bnd_diffusion) then
@@ -569,6 +579,12 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
           enddo ; enddo ; endif
           do j=js,je ; do i=is,ie
             Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j)
+            if (associated(Reg%Tr(m)%diffusionc_xy)) then
+                Reg%Tr(m)%diffusionc_xy(i,j,k) = dTr(i,j) * Idt
+            endif
+            if (associated(Reg%Tr(m)%diffusion_xy)) then
+                Reg%Tr(m)%diffusion_xy(i,j,k) = dTr(i,j) * Idt * (h(i,j,k)+h_neglect)
+            endif
           enddo ; enddo
         enddo
 

src/tracer/MOM_tracer_registry.F90

@@ -53,6 +53,8 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit
                            cmor_name, cmor_units, cmor_longname, net_surfflux_name, NLT_budget_name, &
                            net_surfflux_longname, tr_desc, OBC_inflow, OBC_in_u, OBC_in_v, ad_x, ad_y, &
                            df_x, df_y, ad_2d_x, ad_2d_y, df_2d_x, df_2d_y, advection_xy, registry_diags, &
+                           advectionc_xy, advectionc_x, advectionc_y, &
+                           diffusionc_xy, diffusion_xy, &
                            conc_scale, flux_nameroot, flux_longname, flux_units, flux_scale, &
                            convergence_units, convergence_scale, cmor_tendprefix, diag_form, &
                            restart_CS, mandatory, underflow_conc, Tr_out)
@@ -100,6 +102,12 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit
 
   real, dimension(:,:,:), optional, pointer     :: advection_xy !< convergence of lateral advective tracer fluxes
                                                                 !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
+  real, dimension(:,:,:), optional, pointer     :: advectionc_xy !< convergence of lateral advection
+  real, dimension(:,:,:), optional, pointer     :: diffusionc_xy !< convergence of lateral diffusion
+  real, dimension(:,:,:), optional, pointer     :: diffusion_xy !< convergence of lateral diffusive tracer fluxes
+  real, dimension(:,:,:), optional, pointer     :: advectionc_x !< lateral advection concentration
+  real, dimension(:,:,:), optional, pointer     :: advectionc_y !< lateral advection concentration
+
   logical,              optional, intent(in)    :: registry_diags !< If present and true, use the registry for
                                                                 !! the diagnostics of this tracer.
   real,                 optional, intent(in)    :: conc_scale   !< A scaling factor used to convert the concentration
@@ -245,6 +253,11 @@ subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, unit
   if (present(df_2d_x)) then ; if (associated(df_2d_x)) Tr%df2d_x => df_2d_x ; endif
 
   if (present(advection_xy)) then ; if (associated(advection_xy)) Tr%advection_xy => advection_xy ; endif
+  if (present(advectionc_xy)) then; if (associated(advectionc_xy)) Tr%advectionc_xy => advectionc_xy ; endif
+  if (present(diffusionc_xy)) then; if (associated(diffusionc_xy)) Tr%diffusionc_xy => diffusionc_xy ; endif
+  if (present(diffusion_xy)) then;  if (associated(diffusion_xy)) Tr%diffusion_xy  => diffusion_xy ;  endif
+  if (present(advectionc_x)) then;  if (associated(advectionc_x)) Tr%advectionc_x  => advectionc_x ;  endif
+  if (present(advectionc_y)) then;  if (associated(advectionc_y)) Tr%advectionc_y  => advectionc_y ;  endif
 
   if (present(restart_CS)) then
     ! Register this tracer to be read from and written to restart files.
@@ -434,16 +447,45 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
     if (Tr%id_hbd_dfx_2d > 0) call safe_alloc_ptr(Tr%hbd_dfx_2d,IsdB,IedB,jsd,jed)
     if (Tr%id_hbd_dfy_2d > 0) call safe_alloc_ptr(Tr%hbd_dfy_2d,isd,ied,JsdB,JedB)
 
+    ! The following diagnostics for generic tracer budgets were
+    ! originally developed by Enhui Lao, Fan Yang, and Mathieu Poupon.
+
     Tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", &
         diag%axesTL, Time, &
-        'Horizontal convergence of residual mean advective fluxes of '//trim(lowercase(flux_longname)), &
-        conv_units, v_extensive=.true., conversion=Tr%conv_scale*US%s_to_T)
+        'Horizontal convergence of residual mean advective fluxes of '//&
+        trim(lowercase(flux_longname)), conv_units, v_extensive=.true., &
+        conversion=Tr%conv_scale*US%s_to_T)
+    Tr%id_difc_xy = register_diag_field('ocean_model',trim(shortnm)//"_diffusionc_xy", &
+        diag%axesTL, Time, "Horizontal convergence of residual mean diffusive fluxes of "//&
+        trim(shortnm)//' concentration', trim(units)//' s-1')
+    Tr%id_dif_xy = register_diag_field('ocean_model',trim(shortnm)//"_diffusion_xy", &
+        diag%axesTL, Time, "Horizontal convergence of residual mean diffusive fluxes of "//trim(shortnm), &
+        conv_units, conversion=Tr%conv_scale*US%s_to_T)
+    Tr%id_advc_xy = register_diag_field('ocean_model',trim(shortnm)//"_advectionc_xy", &
+        diag%axesTL, Time, "Horizontal convergence of residual mean advective fluxes of "//&
+        trim(shortnm)//' concentration', trim(units)//' s-1')
+    Tr%id_advc_x  = register_diag_field("ocean_model",trim(shortnm)//'_advectionc_x', &
+        diag%axesTL, Time, "Horizontal x mean advective fluxes of "//trim(shortnm)//' concentration', &
+        trim(units)//' s-1')
+    Tr%id_advc_y  = register_diag_field("ocean_model",trim(shortnm)//'_advectionc_y', &
+        diag%axesTL, Time, "Horizontal y mean advective fluxes of "//trim(shortnm)//' concentration', &
+        trim(units)//' s-1')
     Tr%id_adv_xy_2d = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy_2d", &
         diag%axesT1, Time, &
         'Vertical sum of horizontal convergence of residual mean advective fluxes of '//&
         trim(lowercase(flux_longname)), conv_units, conversion=Tr%conv_scale*US%s_to_T)
     if ((Tr%id_adv_xy > 0) .or. (Tr%id_adv_xy_2d > 0)) &
       call safe_alloc_ptr(Tr%advection_xy,isd,ied,jsd,jed,nz)
+    if (Tr%id_difc_xy > 0) &
+      call safe_alloc_ptr(Tr%diffusionc_xy,isd,ied,jsd,jed,nz)
+    if (Tr%id_dif_xy > 0) &
+      call safe_alloc_ptr(Tr%diffusion_xy,isd,ied,jsd,jed,nz)
+    if (Tr%id_advc_xy > 0) &
+      call safe_alloc_ptr(Tr%advectionc_xy,isd,ied,jsd,jed,nz)
+    if (Tr%id_advc_x > 0) &
+      call safe_alloc_ptr(Tr%advectionc_x,isd,ied,jsd,jed,nz)
+    if (Tr%id_advc_y > 0) &
+      call safe_alloc_ptr(Tr%advectionc_y,isd,ied,jsd,jed,nz)
 
     Tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', &
         diag%axesTL, Time, &
@@ -739,6 +781,11 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag)
     if (Tr%id_dfx_2d > 0) call post_data(Tr%id_dfx_2d, Tr%df2d_x, diag)
     if (Tr%id_dfy_2d > 0) call post_data(Tr%id_dfy_2d, Tr%df2d_y, diag)
     if (Tr%id_adv_xy > 0) call post_data(Tr%id_adv_xy, Tr%advection_xy, diag, alt_h=h_diag)
+    if (Tr%id_advc_xy > 0) call post_data(Tr%id_advc_xy, Tr%advectionc_xy, diag, alt_h=h_diag)
+    if (Tr%id_difc_xy > 0) call post_data(Tr%id_difc_xy, Tr%diffusionc_xy, diag, alt_h=h_diag)
+    if (Tr%id_dif_xy > 0) call post_data(Tr%id_dif_xy, Tr%diffusion_xy, diag, alt_h=h_diag)
+    if (Tr%id_advc_x > 0)  call post_data(Tr%id_advc_x, Tr%advectionc_x, diag, alt_h=h_diag)
+    if (Tr%id_advc_y > 0)  call post_data(Tr%id_advc_y, Tr%advectionc_y, diag, alt_h=h_diag)
     if (Tr%id_adv_xy_2d > 0) then
       work2d(:,:) = 0.0
       do k=1,nz ; do j=js,je ; do i=is,ie

src/tracer/MOM_tracer_types.F90

@@ -48,6 +48,12 @@ module MOM_tracer_types
 
   real, dimension(:,:,:), pointer :: advection_xy   => NULL() !< convergence of lateral advective tracer fluxes
                                                               !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
+  real, dimension(:,:,:), pointer :: advectionc_xy  => NULL() !< convergence of lateral advection concentration
+  real, dimension(:,:,:), pointer :: diffusionc_xy  => NULL() !< convergence of lateral diffusion concentration
+  real, dimension(:,:,:), pointer :: diffusion_xy   => NULL() !< convergence of lateral diffusion content
+  real, dimension(:,:,:), pointer :: advectionc_x   => NULL() !< lateral advection concentration
+  real, dimension(:,:,:), pointer :: advectionc_y   => NULL() !< lateral advection concentration
+
 !  real, dimension(:,:,:), pointer :: diff_cont_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
 !                                                              !! [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1]
 !  real, dimension(:,:,:), pointer :: diff_conc_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
@@ -108,6 +114,8 @@ module MOM_tracer_types
   integer :: id_hbd_dfx_2d = -1, id_hbd_dfy_2d = -1
   integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1
   integer :: id_adv_xy = -1, id_adv_xy_2d = -1
+  integer :: id_advc_xy = -1, id_advc_x = -1, id_advc_y = -1
+  integer :: id_difc_xy = -1, id_dif_xy = -1
   integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1
   integer :: id_hbdxy_cont = -1, id_hbdxy_cont_2d = -1, id_hbdxy_conc = -1
   integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1