Sync to dev/gfdl and re-Added generic tracer budget diagnostics

config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90

@@ -210,7 +210,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, &
@@ -612,13 +612,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

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_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

@@ -401,6 +401,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
@@ -570,6 +580,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_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