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 [H ~> m or kg m-2]
[in]	evap_cfl_limit	Limit on the fraction of the water that can be fluxed out of the top layer in a timestep [nondim]
[in]	minimum_forcing_depth	The smallest depth over which fluxes can be applied [m]
[in]	in_flux_optional	The total time-integrated amount of tracer that enters with freshwater
[in]	out_flux_optional	The total time-integrated amount of tracer that leaves with freshwater
[in]	update_h_opt	Optional flag to determine whether h should be updated

Definition at line 225 of file MOM_tracer_diabatic.F90.

 
  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 [H ~> m or kg m-2]
  real,                                       intent(in   ) :: evap_CFL_limit !< Limit on the fraction of the
                                                                  !! water that can be fluxed out of the top
                                                                  !! layer in a timestep [nondim]
  real,                                       intent(in   ) :: minimum_forcing_depth !< The smallest depth over
                                                                  !! which fluxes can be applied [m]
  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
  logical,                          optional, intent(in) :: update_h_opt  !< Optional flag to determine whether
                                                                  !! h should be updated
 
  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 [Pa s].
  real, dimension(SZI_(G)) :: &
    netMassInOut, &  ! surface water fluxes [H ~> m or kg m-2] over time step
    netMassIn,    &  ! mass entering ocean surface [H ~> m or kg m-2] over a time step
    netMassOut       ! mass leaving ocean surface [H ~> m or kg m-2] 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
 
  ! If no freshwater fluxes, nothing needs to be done in this routine
  if ( (.not. associated(fluxes%netMassIn)) .or. (.not. associated(fluxes%netMassOut)) ) 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 ~> m or kg m-2] 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_diabatic.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
 

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 [H ~> m or kg m-2]
[in]	ea	amount of fluid entrained from the layer above [H ~> m or kg m-2]
[in]	eb	amount of fluid entrained from the layer below [H ~> m or kg m-2]
[in,out]	tr	tracer concentration in concentration units [CU]
[in]	dt	amount of time covered by this call [s]
[in]	sfc_flux	surface flux of the tracer [CU kg m-2 s-1]
[in]	btm_flux	The (negative upward) bottom flux of the tracer [CU kg m-2 s-1]
[in,out]	btm_reservoir	amount of tracer in a bottom reservoir [CU kg m-2]; formerly [CU m]
[in]	sink_rate	rate at which the tracer sinks [m s-1]
[in]	convert_flux_in	True if the specified sfc_flux needs to be integrated in time

Definition at line 27 of file MOM_tracer_diabatic.F90.

  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
                                                                     !! [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: ea     !< amount of fluid entrained from the layer
                                                                     !! above [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: eb     !< amount of fluid entrained from the layer
                                                                     !! below [H ~> m or kg m-2]
  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 [s]
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(in)    :: sfc_flux !< surface flux of the tracer [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 [CU kg m-2 s-1]
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir
                                                                     !! [CU kg m-2]; formerly [CU m]
  real,                             optional,intent(in)    :: sink_rate !< rate at which the tracer sinks [m s-1]
  logical,                          optional,intent(in)    :: convert_flux_in !< True if the specified sfc_flux needs
                                                                     !! to be integrated in time
 
  ! local variables
  real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    sfc_src, &      !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2].
    btm_src         !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2].
  real, dimension(SZI_(G)) :: &
    b1, &           !< b1 is used by the tridiagonal solver [H-1 ~> 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 [nondim].
  real :: h_minus_dsink(SZI_(G),SZK_(GV)) !< The layer thickness minus the
                    !! difference in sinking rates across the layer [H ~> 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 [H ~> m or kg m-2].
  real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2].
  real :: h_tr      !< h_tr is h at tracer points with a h_neglect added to
                    !! ensure positive definiteness [H ~> m or kg m-2].
  real :: h_neglect !< A thickness that is so small it is usually lost
                    !! in roundoff and can be neglected [H ~> m or kg m-2].
  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 (nz == 1) then
    call mom_error(warning, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//&
                            "with only one vertical level")
    return
  endif
 
  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