mom_tracer_diabatic Module Reference

Detailed Description

This module contains routines that implement physical fluxes of tracers (e.g. due to surface fluxes or mixing). These are intended to be called from call_tracer_column_fns in the MOM_tracer_flow_control module.

Functions/Subroutines
subroutine, public	tracer_vertdiff (h_old, ea, eb, dt, tr, G, GV, sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
	This subroutine solves a tridiagonal equation for the final tracer concentrations after the dual-entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking is implemented with an fully implicit upwind advection scheme. More...
subroutine, public	applytracerboundaryfluxesinout (G, GV, Tr, dt, fluxes, h, evap_CFL_limit, minimum_forcing_depth, in_flux_optional, out_flux_optional, update_h_opt)
	This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90 NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface flux of the tracer does not get applied again during a subsequent call to tracer_vertdif. More...

Function/Subroutine Documentation

applytracerboundaryfluxesinout()

subroutine, public mom_tracer_diabatic::applytracerboundaryfluxesinout	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	Tr,
		real, intent(in)	dt,
		type(forcing), intent(in)	fluxes,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h,
		real, intent(in)	evap_CFL_limit,
		real, intent(in)	minimum_forcing_depth,
		real, dimension(szi_(g),szj_(g)), intent(in), optional	in_flux_optional,
		real, dimension(szi_(g),szj_(g)), intent(in), optional	out_flux_optional,
		logical, intent(in), optional	update_h_opt
	)

This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90 NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface flux of the tracer does not get applied again during a subsequent call to tracer_vertdif.

Parameters

[in]	g	Grid structure
[in]	gv	ocean vertical grid structure
[in,out]	tr	Tracer concentration on T-cell
[in]	dt	Time-step over which forcing is applied (s)
[in]	fluxes	Surface fluxes container
[in,out]	h	Layer thickness in H units
[in]	out_flux_optional	Optional flag to determine whether h should be updated

Definition at line 214 of file MOM_tracer_diabatic.F90.

Referenced by advection_test_tracer::advection_test_tracer_column_physics(), boundary_impulse_tracer::boundary_impulse_tracer_column_physics(), ideal_age_example::ideal_age_tracer_column_physics(), isomip_tracer::isomip_tracer_column_physics(), mom_ocmip2_cfc::ocmip2_cfc_column_physics(), and pseudo_salt_tracer::pseudo_salt_tracer_column_physics().

  type(ocean_grid_type),                      intent(in   ) :: g  !< Grid structure
  type(verticalgrid_type),                    intent(in   ) :: gv !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: tr !< Tracer concentration on T-cell
  real,                                       intent(in   ) :: dt !< Time-step over which forcing is applied (s)
  type(forcing),                              intent(in   ) :: fluxes !< Surface fluxes container
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: h  !< Layer thickness in H units
  real,                                       intent(in   ) :: evap_cfl_limit
  real,                                       intent(in   ) :: minimum_forcing_depth
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in   ) :: in_flux_optional ! The total time-integrated amount of tracer!
                                                                             ! that enters with freshwater
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: out_flux_optional ! The total time-integrated amount of tracer!
                                                                              ! that leaves with freshwater
  !< Optional flag to determine whether h should be updated
  logical,                          optional, intent(in) :: update_h_opt

  integer, parameter :: maxgroundings = 5
  integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
  real :: h_limit_fluxes, iforcingdepthscale, idt
  real :: dthickness, dtracer
  real :: fractionofforcing, hold, ithickness
  real :: rivermixconst  ! A constant used in implementing river mixing, in Pa s.
  real, dimension(SZI_(G)) :: &
    netmassinout, &  ! surface water fluxes (H units) over time step
    netmassin,    &  ! mass entering ocean surface (H units) over a time step
    netmassout       ! mass leaving ocean surface (H units) over a time step

  real, dimension(SZI_(G), SZK_(G))                     :: h2d, tr2d
  real, dimension(SZI_(G),SZJ_(G))                      :: in_flux  ! The total time-integrated amount of tracer!
                                                                       ! that enters with freshwater
  real, dimension(SZI_(G),SZJ_(G))                      :: out_flux ! The total time-integrated amount of tracer!
                                                                        ! that leaves with freshwater
  real, dimension(SZI_(G))                              :: in_flux_1d, out_flux_1d
  real                                                  :: hgrounding(maxgroundings)
  real    :: tr_in
  logical :: update_h
  integer :: i, j, is, ie, js, je, k, nz, n, nsw
  character(len=45) :: mesg

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

!  ! Only apply forcing if fluxes%sw is associated.
!  if (.not.ASSOCIATED(fluxes%sw)) return

  in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0
  if(present(in_flux_optional)) then
    do j=js,je ; do i=is,ie
      in_flux(i,j) = in_flux_optional(i,j)
    enddo; enddo
  endif
  if(present(out_flux_optional)) then
    do j=js,je ; do i=is,ie
      out_flux(i,j) = out_flux_optional(i,j)
    enddo ; enddo
  endif

  if (present(update_h_opt)) then
    update_h = update_h_opt
  else
    update_h = .true.
  endif

  idt = 1.0/dt
  numberofgroundings = 0

!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt,    &
!$OMP                                  IforcingDepthScale,minimum_forcing_depth, &
!$OMP                                  numberOfGroundings,iGround,jGround,update_h, &
!$OMP                                  in_flux,out_flux,hGrounding,Idt,evap_CFL_limit) &
!$OMP                          private(h2d,Tr2d,netMassInOut,netMassOut,      &
!$OMP                                  in_flux_1d,out_flux_1d,fractionOfForcing,     &
!$OMP                                  dThickness,dTracer,hOld,Ithickness,           &
!$OMP                                  netMassIn, Tr_in)

  ! Work in vertical slices for efficiency
  do j=js,je

    ! Copy state into 2D-slice arrays
    do k=1,nz ; do i=is,ie
      h2d(i,k) = h(i,j,k)
      tr2d(i,k) = tr(i,j,k)
    enddo ; enddo

    do i = is,ie
      in_flux_1d(i) = in_flux(i,j)
      out_flux_1d(i) = out_flux(i,j)
    enddo
    ! The surface forcing is contained in the fluxes type.
    ! We aggregate the thermodynamic forcing for a time step into the following:
    ! These should have been set and stored during a call to applyBoundaryFluxesInOut
    ! netMassIn    = net mass entering at ocean surface over a timestep
    ! netMassOut   = net mass leaving ocean surface (H units) over a time step.
    !                netMassOut < 0 means mass leaves ocean.

    ! Note here that the aggregateFW flag has already been taken care of in the call to
    ! applyBoundaryFluxesInOut
    do i=is,ie
        netmassout(i) = fluxes%netMassOut(i,j)
        netmassin(i)  = fluxes%netMassIn(i,j)
    enddo

    ! Apply the surface boundary fluxes in three steps:
    ! A/ update concentration from mass entering the ocean
    ! B/ update concentration from mass leaving ocean.
    do i=is,ie
      if (g%mask2dT(i,j)>0.) then

        ! A/ Update tracer due to incoming mass flux.
        do k=1,1

          ! Change in state due to forcing
          dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
          dtracer = 0.

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
          netmassin(i) = netmassin(i) - dthickness
          dtracer = dtracer + in_flux_1d(i)
          in_flux_1d(i) = 0.0

          ! Update state
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.0) then
            ithickness  = 1.0/h2d(i,k)      ! Inverse new thickness
            ! The "if"s below avoid changing T/S by roundoff unnecessarily
            if (dthickness /= 0. .or. dtracer /= 0.) tr2d(i,k) = (hold*tr2d(i,k)+ dtracer)*ithickness
          endif

        enddo ! k=1,1

        ! B/ Update tracer from mass leaving ocean
        do k=1,nz

          ! Place forcing into this layer if this layer has nontrivial thickness.
          ! For layers thin relative to 1/IforcingDepthScale, then distribute
          ! forcing into deeper layers.
          iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth*gv%m_to_H - netmassout(i) )
          ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
          fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)

          ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
          ! limit the forcing applied to this cell, leaving the remaining forcing to
          ! be distributed downwards.
          if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
            fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
          endif

          ! Change in state due to forcing
          dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
          ! Note this is slightly different to how salt is currently treated
          dtracer =  fractionofforcing*out_flux_1d(i)

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then new netMassOut vanishes.
          netmassout(i) = netmassout(i) - dthickness
          out_flux_1d(i) = out_flux_1d(i) - dtracer

          ! Update state by the appropriate increment.
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.) then
            ithickness  = 1.0/h2d(i,k) ! Inverse of new thickness
            tr2d(i,k)    = (hold*tr2d(i,k) + dtracer)*ithickness
          endif

        enddo ! k

      endif

      ! If anything remains after the k-loop, then we have grounded out, which is a problem.
      if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0) then
!$OMP critical
        numberofgroundings = numberofgroundings +1
        if (numberofgroundings<=maxgroundings) then
          iground(numberofgroundings) = i ! Record i,j location of event for
          jground(numberofgroundings) = j ! warning message
          hgrounding(numberofgroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i))
        endif
!$OMP end critical
      endif

    enddo ! i

    ! Step C/ copy updated tracer concentration from the 2d slice now back into model state.
    do k=1,nz ; do i=is,ie
      tr(i,j,k) = tr2d(i,k)
    enddo ; enddo

    if (update_h) then
      do k=1,nz ; do i=is,ie
        h(i,j,k) = h2d(i,k)
      enddo ; enddo
    endif

  enddo ! j-loop finish

  if (numberofgroundings>0) then
    do i = 1, min(numberofgroundings, maxgroundings)
      write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
                             g%geoLatT( iground(i), jground(i)) , hgrounding(i)
      call mom_error(warning, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
                              "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.)
    enddo

    if (numberofgroundings - maxgroundings > 0) then
      write(mesg, '(i4)') numberofgroundings - maxgroundings
      call mom_error(warning, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
                              trim(mesg) // " groundings remaining")
    endif
  endif

Here is the caller graph for this function:

tracer_vertdiff()

subroutine, public mom_tracer_diabatic::tracer_vertdiff	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	h_old,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	ea,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	eb,
		real, intent(in)	dt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)	tr,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional	sfc_flux,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional	btm_flux,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout), optional	btm_reservoir,
		real, intent(in), optional	sink_rate,
		logical, intent(in), optional	convert_flux_in
	)

This subroutine solves a tridiagonal equation for the final tracer concentrations after the dual-entrainments, and possibly sinking or surface and bottom sources, are applied. The sinking is implemented with an fully implicit upwind advection scheme.

Parameters

[in]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h_old	layer thickness before entrainment (m or kg m-2)
[in]	ea	amount of fluid entrained from the layer above (units of h_work)
[in]	eb	amount of fluid entrained from the layer below (units of h_work)
[in,out]	tr	tracer concentration (in concentration units CU)
[in]	dt	amount of time covered by this call (seconds)
[in]	sfc_flux	surface flux of the tracer (in CU * kg m-2 s-1)
[in]	btm_flux	The (negative upward) bottom flux of the tracer, in units of (CU * kg m-2 s-1)
[in,out]	btm_reservoir	amount of tracer in a bottom reservoir (units of CU kg m-2; formerly CU m)
[in]	sink_rate	rate at which the tracer sinks, in m s-1
[in]	convert_flux_in	True if the specified sfc_flux needs to be integrated in time

Definition at line 26 of file MOM_tracer_diabatic.F90.

Referenced by advection_test_tracer::advection_test_tracer_column_physics(), boundary_impulse_tracer::boundary_impulse_tracer_column_physics(), mom_diabatic_driver::diabatic(), ideal_age_example::ideal_age_tracer_column_physics(), isomip_tracer::isomip_tracer_column_physics(), mom_ocmip2_cfc::ocmip2_cfc_column_physics(), and pseudo_salt_tracer::pseudo_salt_tracer_column_physics().

  type(ocean_grid_type),                     intent(in)    :: g             !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv            !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_old         !< layer thickness before entrainment (m or kg m-2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: ea            !< amount of fluid entrained from the layer above (units of h_work)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: eb            !< amount of fluid entrained from the layer below (units of h_work)
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr            !< tracer concentration (in concentration units CU)
  real,                                      intent(in)    :: dt            !< amount of time covered by this call (seconds)
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(in)    :: sfc_flux      !< surface flux of the tracer (in CU * kg m-2 s-1)
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(in)    :: btm_flux      !< The (negative upward) bottom flux of the tracer,
                                                                            !! in units of (CU * kg m-2 s-1)
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir (units of CU kg m-2; formerly CU m)
  real,                             optional,intent(in)    :: sink_rate     !< rate at which the tracer sinks, in m s-1
  logical,                          optional,intent(in)  :: convert_flux_in    !< True if the specified sfc_flux needs to be integrated in time

  real :: sink_dist ! The distance the tracer sinks in a time step, in m or kg m-2.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    sfc_src, &      ! The time-integrated surface source of the tracer, in
                    ! units of m or kg m-2 times a concentration.
    btm_src         ! The time-integrated bottom source of the tracer, in
                    ! units of m or kg m-2  times a concentration.
  real, dimension(SZI_(G)) :: &
    b1, &           ! b1 is used by the tridiagonal solver, in m-1 or m2 kg-1.
    d1              ! d1=1-c1 is used by the tridiagonal solver, nondimensional.
  real :: c1(szi_(g),szk_(gv))    ! c1 is used by the tridiagonal solver, ND.
  real :: h_minus_dsink(szi_(g),szk_(gv)) ! The layer thickness minus the
                    ! difference in sinking rates across the layer, in m or kg m-2.
                    ! By construction, 0 <= h_minus_dsink < h_work.
  real :: sink(szi_(g),szk_(gv)+1) ! The tracer's sinking distances at the
                    ! interfaces, limited to prevent characteristics from
                    ! crossing within a single timestep, in m or kg m-2.
  real :: b_denom_1 ! The first term in the denominator of b1, in m or kg m-2.
  real :: h_tr      ! h_tr is h at tracer points with a h_neglect added to
                    ! ensure positive definiteness, in m or kg m-2.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in m.
  logical :: convert_flux = .true.


  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  if (present(convert_flux_in)) convert_flux = convert_flux_in
  h_neglect = gv%H_subroundoff
  sink_dist = 0.0
  if (present(sink_rate)) sink_dist = (dt*sink_rate) * gv%m_to_H
!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, &
!$OMP                               sink_rate,btm_reservoir,nz,sink_dist,ea,      &
!$OMP                               h_old,convert_flux,h_neglect,eb,tr) &
!$OMP                       private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
!$OMP do
  do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo; enddo
  if (present(sfc_flux)) then
    if(convert_flux) then
!$OMP do
      do j = js, je; do i = is,ie
        sfc_src(i,j) = (sfc_flux(i,j)*dt) * gv%kg_m2_to_H
      enddo; enddo
    else
!$OMP do
      do j = js, je; do i = is,ie
        sfc_src(i,j) = sfc_flux(i,j)
      enddo; enddo
    endif
  endif
  if (present(btm_flux)) then
    if(convert_flux) then
!$OMP do
      do j = js, je; do i = is,ie
        btm_src(i,j) = (btm_flux(i,j)*dt) * gv%kg_m2_to_H
      enddo; enddo
    else
!$OMP do
      do j = js, je; do i = is,ie
        btm_src(i,j) = btm_flux(i,j)
      enddo; enddo
    endif
  endif

  if (present(sink_rate)) then
!$OMP do
    do j=js,je
      ! Find the sinking rates at all interfaces, limiting them if necesary
      ! so that the characteristics do not cross within a timestep.
      !   If a non-constant sinking rate were used, that would be incorprated
      ! here.
      if (present(btm_reservoir)) then
        do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
        do k=2,nz ; do i=is,ie
          sink(i,k) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
        enddo ; enddo
      else
        do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
        ! Find the limited sinking distance at the interfaces.
        do k=nz,2,-1 ; do i=is,ie
          if (sink(i,k+1) >= sink_dist) then
            sink(i,k) = sink_dist
            h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,k+1) - sink(i,k))
          elseif (sink(i,k+1) + h_old(i,j,k) < sink_dist) then
            sink(i,k) = sink(i,k+1) + h_old(i,j,k)
            h_minus_dsink(i,k) = 0.0
          else
            sink(i,k) = sink_dist
            h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,k+1)) - sink(i,k)
          endif
        enddo ; enddo
      endif
      do i=is,ie
        sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
      enddo

      ! Now solve the tridiagonal equation for the tracer concentrations.
      do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = b_denom_1 * b1(i)
        h_tr = h_old(i,j,1) + h_neglect
        tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
      endif ; enddo
      do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        c1(i,k) = eb(i,j,k-1) * b1(i)
        b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,k)) + &
                    h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        h_tr = h_old(i,j,k) + h_neglect
        tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
                             (ea(i,j,k) + sink(i,k)) * tr(i,j,k-1))
      endif ; enddo ; enddo
      do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        c1(i,nz) = eb(i,j,nz-1) * b1(i)
        b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
                    h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
        h_tr = h_old(i,j,nz) + h_neglect
        tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                              (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
      endif ; enddo
      if (present(btm_reservoir)) then ; do i=is,ie ; if (g%mask2dT(i,j)>0.5) then
        btm_reservoir(i,j) = btm_reservoir(i,j) + &
                             (sink(i,nz+1)*tr(i,j,nz)) * gv%H_to_kg_m2
      endif ; enddo ; endif

      do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
      endif ; enddo ; enddo
    enddo
  else
!$OMP do
    do j=js,je
      do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        h_tr = h_old(i,j,1) + h_neglect
        b_denom_1 = h_tr + ea(i,j,1)
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = h_tr * b1(i)
        tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
       endif
      enddo
      do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        c1(i,k) = eb(i,j,k-1) * b1(i)
        h_tr = h_old(i,j,k) + h_neglect
        b_denom_1 = h_tr + d1(i) * ea(i,j,k)
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1))
      endif ; enddo ; enddo
      do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        c1(i,nz) = eb(i,j,nz-1) * b1(i)
        h_tr = h_old(i,j,nz) + h_neglect
        b_denom_1 = h_tr + d1(i)*ea(i,j,nz)
        b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz))
        tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                              ea(i,j,nz) * tr(i,j,nz-1))
      endif ; enddo
      do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
      endif ; enddo ; enddo
    enddo
  endif

!$OMP end parallel

Here is the caller graph for this function: