Apply changes for BGC budgets

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

src/tracer/MOM_generic_tracer.F90

@@ -601,9 +601,23 @@ 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)
+            !liao
+            g_tracer%boundary_forcing_tend(i,j,k) = 0
+            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
+            !liao
           enddo ; enddo ; enddo
           call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, &
                             fluxes, h_work, evap_CFL_limit, minimum_forcing_depth)
+            !liao
+            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
+            !liao
         endif
 
          !traverse the linked list till hit NULL

src/tracer/MOM_tracer_advect.F90

@@ -187,9 +187,50 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first
   ! initialize diagnostic fluxes and tendencies
   !$OMP do
   do m=1,ntr
+!    if (associated(Tr(m)%ad_x)) then
+!      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+!        Tr(m)%ad_x(I,j,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    if (associated(Tr(m)%ad_y)) then
+!      do k=1,nz ; do J=jsd,jed ; do i=isd,ied
+!        Tr(m)%ad_y(i,J,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    if (associated(Tr(m)%advection_xy)) then
+!      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+!        Tr(m)%advection_xy(i,j,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    !liao
+!    if (associated(Tr(m)%advectionc_xy)) then
+!      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+!        Tr(m)%advectionc_xy(i,j,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    if (associated(Tr(m)%advectionc_x)) then
+!      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+!        Tr(m)%advectionc_x(i,j,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    if (associated(Tr(m)%advectionc_y)) then
+!      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+!        Tr(m)%advectionc_y(i,j,k) = 0.0
+!      enddo ; enddo ; enddo
+!    endif
+!    !liao
+!    if (associated(Tr(m)%ad2d_x)) then
+!      do j=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_x(I,j) = 0.0 ; enddo ; enddo
+!    endif
+!    if (associated(Tr(m)%ad2d_y)) then
+!      do J=jsd,jed ; do i=isd,ied ; Tr(m)%ad2d_y(i,J) = 0.0 ; enddo ; enddo
+!    endif
     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
@@ -362,6 +403,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 !liao
   real, dimension(SZIB_(G)) :: &
     hlst, &             ! Work variable [H L2 ~> m3 or kg].
     Ihnew, &            ! Work variable [H-1 L-2 ~> m-3 or kg-1].
@@ -652,6 +694,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)  !liao
             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
@@ -671,6 +714,18 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
                                           Idt * G%IareaT(i,j)
         endif ; enddo
       endif
+      !liao
+      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
+      !liao
 
     enddo
 
@@ -733,6 +788,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  !liao
   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
@@ -1039,6 +1095,7 @@ 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) !liao
         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
@@ -1051,6 +1108,18 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
                                           G%IareaT(i,j)
         endif ; enddo
       endif
+      !liao
+      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
+      !liao
 
     enddo
   endif ; enddo ! End of j-loop.

src/tracer/MOM_tracer_hor_diff.F90

@@ -400,6 +400,13 @@ 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
+    ! mpoupon
+    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
+    ! mpoupon
   enddo
 
   if (CS%use_hor_bnd_diffusion) then
@@ -569,6 +576,11 @@ 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)
+            ! mpoupon
+            if (associated(Reg%Tr(m)%diffusionc_xy)) then
+                Reg%Tr(m)%diffusionc_xy(i,j,k) = dTr(i,j) * Idt
+            endif
+            ! mpoupon
           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, & !liao
+                           diffusionc_xy, & ! mpoupon
                            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,11 @@ 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 !liao
+  real, dimension(:,:,:), optional, pointer     :: diffusionc_xy !< convergence of lateral diffusion !mpoupon
+  real, dimension(:,:,:), optional, pointer     :: advectionc_x !< lateral advection concentration !liao
+  real, dimension(:,:,:), optional, pointer     :: advectionc_y !< lateral advection concentration !liao
+
   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 +252,10 @@ 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 !liao  
+  if (present(diffusionc_xy)) then; if (associated(diffusionc_xy)) Tr%diffusionc_xy => diffusionc_xy ; endif !mpoupon
+  if (present(advectionc_x)) then;  if (associated(advectionc_x)) Tr%advectionc_x  => advectionc_x ;  endif !liao
+  if (present(advectionc_y)) then;  if (associated(advectionc_y)) Tr%advectionc_y  => advectionc_y ;  endif !liao
 
   if (present(restart_CS)) then
     ! Register this tracer to be read from and written to restart files.
@@ -436,14 +447,43 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u
 
     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)
+    ! mpoupon
+    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), &
+        trim(units)//' s-1')
+    ! mpoupon
+    !liao
+    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), &
+        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), &
+        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), &
+        trim(units)//' s-1')
+    !liao
     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)
+    !mpoupon
+    if (Tr%id_difc_xy > 0) &
+      call safe_alloc_ptr(Tr%diffusionc_xy,isd,ied,jsd,jed,nz)
+    !mpoupon
+    !liao
+    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)
+    !liao
 
     Tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', &
         diag%axesTL, Time, &
@@ -739,6 +779,10 @@ 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) !liao
+    if (Tr%id_difc_xy > 0) call post_data(Tr%id_difc_xy, Tr%diffusionc_xy, diag, alt_h=h_diag) !mpoupon
+    if (Tr%id_advc_x > 0)  call post_data(Tr%id_advc_x, Tr%advectionc_x, diag, alt_h=h_diag) !liao
+    if (Tr%id_advc_y > 0)  call post_data(Tr%id_advc_y, Tr%advectionc_y, diag, alt_h=h_diag) !liao
     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,11 @@ 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 liao
+  real, dimension(:,:,:), pointer :: diffusionc_xy  => NULL() !< convergence of lateral diffusion concentration mpoupon
+  real, dimension(:,:,:), pointer :: advectionc_x  => NULL()  !< lateral advection concentration liao
+  real, dimension(:,:,:), pointer :: advectionc_y  => NULL()  !< lateral advection concentration liao
+
 !  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 +113,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 !liao
+  integer :: id_difc_xy = -1 ! mpoupon
   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