tracer diagnostics

Message

nbanas · #1 Unread post by **nbanas** » Tue Nov 08, 2005 8:10 pm

hi,
I'm trying to understand the relationship between the five diagnostics one gets for each tracer when DIAGNOSTICS_TS is turned on. I would have guessed that (for example)

salt_rate = salt_hdiff + salt_vdiff + salt_hadv + salt_vadv

(or maybe that some of these terms were on the other side of the equal sign), and that since this is a linear budget, that this would be strictly true for any averaging period, and cell by cell. It doesn't seem to be, though: in many places these don't balance cell by cell within a factor of 2. Could it have to do with correlations with layer thickness over the tidal cycle (I'm getting diagnostics every 25 h)?

Also--for biological tracers, does tracer_rate include the net change wrought by fasham.h, or is it only the physics?

thanks,
Neil

wilkin · #2 Unread post by **wilkin** » Mon Nov 14, 2005 2:39 pm

Your question is timely: Motivated by some of the group discussion at the User's meeting last month, I had prepared a few snipets of documentation for the tracer diagnostics to go into the web pages, but to answer your question I'll post them here:

In short, the balance you propose is correct. I just double-checked one of my CBLAST model runs and get closure of temp/salt_rate = hadv+vadv+vdiff (I had no hdiff term) to within 0.01%. Either the error is in the saved hdiff term, or somewhere else in your configuration.

I believe the net effects of the internal transformations in the biological code will end up in the _rate term. If you look in the code you'lll see the term is computed simply as the change from the preceeding time step.

John.

Documentation on tracer diagnostic conventions:

The sign convention for the diagnostics terms is that all terms have been taken over to the right-hand-side, *except* for the time rate of change itself.

Therefore, a check on the sum of any set of tracer diagnostics should find that:

Tname_rate = Tname_hadv + Tname_vadv + Tname_hdiff + Tname_vdiff

To reconcile the the vertical integral of the vertical diffusion term with net surface heat flux (variable 'shflux' in the output history and averages file), consider that:

integral_-h^zeta d/dz(kappa_t dT/dz) dz = stflx(surface)-stflx(-h)

where stflx are the vertical boundary conditions on vertical diffusion, i.e. the surface and bottom *heat* fluxes Q(surface) and Q(-h) expressed in units of degC m s^-1.

Assuming Q(-h)=0, then the vertical integral of the vertical diffusion term will be positive if the net air-sea flux is warming the ocean, and negative if it is cooling the ocean.

The netcdf attributes for shflux state that:
shflux:negative_value = "upward flux, cooling"
shflux:positive_value = "downward flux, warming"

so you should find that these terms have the same sign.

If you divide shflux (W m^-2) by rho0*cp (1025*3985 degC m3 s^-1) to get the same units as the diagnostics (degC m s^-1) you should obtain the same numeric value as the vertical integral of temp_vdiff (to within the accuracy of the vertical integration step, which requires you make some assumption about the average zeta over the diagostics averaging interval when computing the vertical layer thicknesses).

wmartin · #3 Unread post by **wmartin** » Sat Feb 04, 2006 3:38 am

I have been having similar problems with the T and S diagnostics. After looking at the code in step3d_t_tile, I think there may be a problem when the MPDATA option is set. In that case it looks like the basic horizontal and vertical advection terms are computed and summed without any calculation of the diagnostics. After that the "anti-difussive" corrections are made in the MPDATA logic and then the diagnostic terms are set equal to the corrections. Thus, when it's all done, the diagnostics represent the corrections rather than the actual values of the underlying terms.

This is extremely complex code at the *.F level with all the cpp options. I tried to sort it out, but gave up. Please let me know if I am on the right track and if so how to fix it.

Thanks
Wayne

wilkin · #4 Unread post by **wilkin** » Sat Feb 04, 2006 5:19 pm

I don't know for sure that anyone has tested the diagnostics terms with the MPDATA option. I certainly haven't.

If you're having trouble following the *.F code, then browse the *.f90 for your configuration and follow whether it's doing the right things. The diagnostics should certainly balance regardless of the advection option, but MPDATA is a new option so it will take some debugging.

When you figure out the problem, please make sure the information gets back to us.

John.

cae · #5 Unread post by **cae** » Mon Feb 13, 2006 6:59 pm

Hi.

I was able to get the tracer budgets to work along with MPDATA by making 4 changes to step3d_t.F. It seemed to me that the diagnostic tracer work arrays were not being fully updated for MPDATA. I'm not 100% sure this is right, but I believe it's consistent. The lines are listed below between comments with my initials (CAE). I've included enough of the code to identify where in step3d_t.F these changes are needed

Code: Select all

    686 #  ifdef TS_MPDATA
    687               Ta(i,j,k,itrc)=t(i,j,k,3,itrc)-cff1
    688 !CAE
    689               DiaTwrk(i,j,k,itrc,iThadv)=-cff1
    690 !END CAE
    691 #  else
    692               t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
    693 #   ifdef DIAGNOSTICS_TS
    694               DiaTwrk(i,j,k,itrc,iThadv)=-cff1
    695 #   endif
    696 #  endif
    697             END DO
    698           END DO
    699         END DO K_LOOP
    700       END DO T_LOOP
    701 !
    702 !-----------------------------------------------------------------------    
    703 !  Time-step vertical advection term.
    704 !---------------------------------------------------------------------

Code: Select all

    848 !
    849 !  Time-step vertical advection term.
    850 # ifdef DIAGNOSTICS_TS
    851 !  Convert units of tracer diagnostic terms to Tunits.
    852 # endif
    853 !
    854           DO i=I_RANGE
    855             CF(i,0)=dt(ng)*pm(i,j)*pn(i,j)
    856           END DO
    857           DO k=1,N(ng)
    858             DO i=I_RANGE
    859               cff1=CF(i,0)*(FC(i,k)-FC(i,k-1))
    860 # ifdef TS_MPDATA
    861               Ta(i,j,k,itrc)=(Ta(i,j,k,itrc)-cff1)*oHz(i,j,k)
    862 !CAE
    863               DiaTwrk(i,j,k,itrc,iTvadv)=-cff1
    864 !END CAE
    865 # else
    866               t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1

Code: Select all

    976 !
    977 !  Time-step corrected horizontal advection (Tunits m).
    978 !
    979           DO j=Jstr,Jend
    980             DO i=Istr,Iend
    981               cff1=dt(ng)*(FX(i+1,j)-FX(i,j)+         &
    982      &                     FE(i,j+1)-FE(i,j))*         &
    983      &                     pm(i,j)*pn(i,j)
    984               t(i,j,k,nnew,itrc)=Ta(i,j,k,itrc)*Hz(i,j,k)-cff1
    985 #  ifdef DIAGNOSTICS_TS
    986 !CAE
    987               DiaTwrk(i,j,k,itrc,iThadv)=DiaTwrk(i,j,k,itrc,iThadv)-         &
    988      &                                   cff1
    989 !END CAE
    990 !ORIGINAL WAY              DiaTwrk(i,j,k,itrc,iThadv)=-cff1
    991 #  endif

Code: Select all

   1011 !
   1012 !  Time-step corrected vertical advection (Tunits).
   1013 #  ifdef DIAGNOSTICS_TS
   1014 !  Convert units of tracer diagnostic terms to Tunits.
   1015 #  endif
   1016 !
   1017           DO i=Istr,Iend
   1018             CF(i,0)=dt(ng)*pm(i,j)*pn(i,j)
   1019           END DO
   1020           DO k=1,N(ng)
   1021             DO i=Istr,Iend
   1022               cff1=CF(i,0)*(FC(i,k)-FC(i,k-1))
   1023               t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
   1024 #  ifdef DIAGNOSTICS_TS
   1025 !CAE
   1026               DiaTwrk(i,j,k,itrc,iTvadv)=DiaTwrk(i,j,k,itrc,iTvadv)-         &
   1027      &                                   cff1
   1028 !END CAE
   1029 !ORIGINAL WAY              DiaTwrk(i,j,k,itrc,iTvadv)=-cff1
   1030               DiaTwrk(i,j,k,itrc,iThadv)=DiaTwrk(i,j,k,itrc,iThadv)*         &
   1031      &                                   oHz(i,j,k)

wmartin · #6 Unread post by **wmartin** » Tue Feb 21, 2006 2:07 am

The code posted by CAE above matches what I came up with and it works for me as well. There are a couple of places in the .F file where it also needs #ifdef DIAGNOSTIS_TS so it won't get included if the diagnostics are not on. In working on this problem I also found that the diagnostics don't reset for NDIA=1. That can be fixed by adding " .or. (nDIA(ng).eq. 1) " to the test at the start of set_diags.F
wayne

Ocean Modeling Discussion