Suspended sediments: nstp or nnew?

[csherwood]

This question demonstrates a fundamental lack of understanding of how ROMS steps tracers.

We are developing a floc model that moves mass from one suspended-sediment tracer class to another. This is a zero-dimensional calculation, performed at every i,j,k location, and integrated over one baroclinic time step.

This happens in our new routine sed_flocmod that is called from sediment.F right before sed_settling

Code: Select all

!***********************************************************************
      SUBROUTINE sediment (ng, tile)
!***********************************************************************
!
      USE sed_bed_mod, ONLY : sed_bed
      USE sed_fluxes_mod, ONLY : sed_fluxes
      USE sed_settling_mod, ONLY : sed_settling
      USE sed_flocs_mod, ONLY : sed_flocmod
      USE sed_surface_mod, ONLY : sed_surface
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!-----------------------------------------------------------------------
!  Compute sediment flocculation
!-----------------------------------------------------------------------
!
      CALL sed_flocmod (ng, tile)
!
!-----------------------------------------------------------------------
!  Compute sediment vertical settling.
!-----------------------------------------------------------------------
!
      CALL sed_settling (ng, tile)
!
!-----------------------------------------------------------------------
!  Compute bed-water column exchanges: erosion and deposition.
!-----------------------------------------------------------------------
!
      CALL sed_fluxes (ng, tile)
!
!-----------------------------------------------------------------------
!  Compute sediment bed stratigraphy.
!-----------------------------------------------------------------------
!
      CALL sed_bed (ng, tile)
!

Why can't we make a local copy of the tracer array, have our way with it, and write it back to the same place? Like this:

Code: Select all

!
!--------------------------------------------------------------------------
!
      J_LOOP : DO j=Jstr,Jend
!
!  Extract mud variables from tracer arrays, place them into
!  local array
      DO ised=1,NCS
         indx = idsed(ised)
         DO k=1,N(ng)
            DO i=Istr,Iend
               susmud(i,k,ised)=t(i,j,k,nstp,indx)
!               susmud(i,k,ised)=t(i,j,k,nnew,indx)?
            ENDDO
         ENDDO
      ENDDO
!-----------------------------------------------------------------------
!  Do zero-d floc calculations here. Move mass from one size class to
!  another, but no net change in mass.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!  Update global tracer variables.
!-----------------------------------------------------------------------
!
      DO ised=1,NCS
         indx = idsed(ised)
         DO k=1,N(ng)
           DO i=Istr,Iend
              t(i,j,k,nstp,indx)=susmud(i,k,ised)
!              t(i,j,k,nnew,indx)=susmud(i,k,ised)
           ENDDO
         ENDDO
      ENDDO

      END DO J_LOOP

When we do this with no actual computations between reading and writing to the tracer arrays, the initial value of one of the tracer classes decays exponentially to zero. Fake values we make up and write to the arrays appear, but are attenuated.

Alternatively, we could read from nstp and write to nnew, like sed_settling does, but that does not appear to work either.

Obviously, I don't know what I am doing with nstp and nnew...how do these work?