Separating horizontal advection term for tracer diagnostics

Message

schen · #1 Unread post by **schen** » Tue Apr 08, 2008 7:58 pm

Hi all,

In ROMS standard TS_DIAGNOSTICS output, the advection terms (x and y) are combined into one term (salt_hadv and temp_adv). I'm trying to separate them in step3d_t.F but am not sure:

1. whether my modifications (see below) are correct or not?

Another related question is that the advection terms are computed in a flux-form. For example, x-advection term is computed as FX (i+1)- FX(i) where FX = u * S .

2. How can I separate d/dx ( u * S ) into two terms ?

Any suggestion is appreciated!!!

-----------------------------------------------------------------------------------------------------
My modifications with TS_MPDATA in step3d_t.F (under flag HADV_SEP)

Code: Select all

!
!  Time-step horizontal advection term.
!
          DO j=J_RANGE
            DO i=I_RANGE
              cff=dt(ng)*pm(i,j)*pn(i,j)
              cff1=cff*(FX(i+1,j)-FX(i,j)+                              &
     &                  FE(i,j+1)-FE(i,j))
#  ifdef HADV_SEP
              cff_x=cff*(FX(i+1,j)-FX(i,j))
              cff_y=cff*(FE(i,j+1)-FE(i,j))
#   endif

#  ifdef TS_MPDATA
              Ta(i,j,k,itrc)=t(i,j,k,3,itrc)-cff1
#  else
              t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
#  endif
#  ifdef DIAGNOSTICS_TS
              DiaTwrk(i,j,k,itrc,iThadv)=-cff1
#  ifdef HADV_SEP
              DiaTwrk(i,j,k,itrc,iTxadv)=-cff_x
              DiaTwrk(i,j,k,itrc,iTyadv)=-cff_y
#   endif
#  endif
            END DO
          END DO

Code: Select all

!
!  Time-step corrected horizontal advection (Tunits m).
!
          DO j=Jstr,Jend
            DO i=Istr,Iend
              cff1=dt(ng)*(FX(i+1,j)-FX(i,j)+                           &
     &                     FE(i,j+1)-FE(i,j))*                          &
     &                     pm(i,j)*pn(i,j)
              t(i,j,k,nnew,itrc)=Ta(i,j,k,itrc)*Hz(i,j,k)-cff1
#  ifdef HADV_SEP
              cff_x=dt(ng)*(FX(i+1,j)-FX(i,j))*                         &
     &                     pm(i,j)*pn(i,j)
              cff_y=dt(ng)*(FE(i,j+1)-FE(i,j))*                         &
     &                     pm(i,j)*pn(i,j)
#  endif

#  ifdef DIAGNOSTICS_TS
              DiaTwrk(i,j,k,itrc,iThadv)=DiaTwrk(i,j,k,itrc,iThadv)-    &
     &                                   cff1

#  ifdef HADV_SEP
              DiaTwrk(i,j,k,itrc,iTxadv)=DiaTwrk(i,j,k,itrc,iTxadv)-    &
     &                                   cff_x
              DiaTwrk(i,j,k,itrc,iTyadv)=DiaTwrk(i,j,k,itrc,iTyadv)-    &
     &                                   cff_y
#  endif

#  endif
            END DO
          END DO

milena · #2 Unread post by **milena** » Wed Apr 09, 2008 1:05 am

Hi schen,

I just did exactly the same thing (not with TS_MPDATA, but it shouldn't make a difference). You don't mention in your post if you modified other routines or not, but there are another few places where you need to tell ROMS that DiaTwrk has two additional diagnostic tracers:
1) in mod_param.F: need to update NDT (number of diagnostic fields)
2) in mod_ncparam.F: in two places you need to tell ROMS about the new variables which will be stored in the netcdf file. Look for iThadv to identify these two places in the routine and include your 2 new cases.
3) in mod_scalars.F: need to define the new diagnostic fields indeces (what you called iTxadv and iTyadv). You need to declare them (similar to "integer :: iThadv") and to define them (look for the place where iThadv is defined).
4) finally, you need to update your varinfo.dat file, including the information for the new variables (look for the case '_hadv' and add two similar cases).
I hope this helps.

I am not sure I understand your second question. You are wondering whether you can get u*S directly instead of the flux?

Milena

schen · #3 Unread post by **schen** » Wed Apr 09, 2008 1:38 am

Hi Milena,
Thanks. I did modify the files you mention

My second question is how to separate d/dx (u x S) into u dS/dx + S du/dx in the code. If you add the advection terms, the Sdu/dx+Sdv/dy+Sdw/dz = 0. But I want to compare u dS/dx and v dS/dy in tracer transport equation. So I'm trying to compute these two separately. Please let me know if you have any suggestion.