[MITgcm-support] Tracer budgets

[Hannah Kleppin]

Dear MITgcm users and developers,

When trying to close elemental budgets in a regional modelling set-up 
(using OBCS and RBCS) using a modified version of the DIC pkg, I came up 
with same general questions regarding the DIC package and Ptracers in 
combination with the OBCS pakage.

1.) Time stepping of the DIC-biogeochemistry model:

 From the documentation and the code I'm not entirely sure what the time 
stepping method of the biogeochemistry is. I would assume that it's 
supposedly a simple Euler forward method. However, I have the impression 
that the calculation of the biogeochemical source and sink terms is done 
with a modified tracer values ( I have difficulties locating it in time; 
I'm using the staggered baroclinic time-stepping in combination with 
advection scheme 77 (which switches of the AB- time stepping for 
thermodynamics, as far as I understand) ), i.e., the tracers are already 
updated for diffusion and advection contributions ( calls to  
timestep_tracer (L420) and cycle_tracer (L498) in ptracers_integrate.F), 
before gchem_forcing_sep.F is called.

In the documentation I found the following sentence regarding 
multi-dimensional advection (section 2.17.2.4)

"In order to incorporate this method into the general model algorithm, 
we compute the effective tendency rather than update the tracer so that 
other terms such as diffusion are using the nthtime-level and not the 
updated n+3/3

quantities".


The full tracer equation should be sth like

Tr^(n+1)=Tr^n+GTr_adv^(n+1/2)*dt+GTr^n_diff*dt +GTr_bio^n *dt

the biological sink and source tendencies (GTr_bio) should then be 
calculated based on TR^n

However, I believe the Tracer when entering subroutine 
dic_biotic_forcing.F is

Tr*:= Tr^n+GTr_adv^(n+1/2)*dt+GTr_diff^n*dt

If I calculate my biogeochemical source and sink terms based on the 
updated (advected and diffused) tracer field, wouldn't there be the need 
to do some kind of time splitting. Or is this done, and I just cannot 
see it?

Or otherwise shouldn't dic_biotic_forcing.F  use both, the unmodified 
Tr^n values and Tr*, with GTr_bio^n= a*Tr^n- b*Tr^n and then integrate 
as Tr^(n+1)= Tr*+GTr_bio^n*dt


2.) OBCS and ptracers

  I calculated the total change of phosphorous in my domain in two 
different ways:

First way: I cumulate the amount of phosphorous in each grid cell of the 
model at certain time intervals and get the change of phosphorus content 
in the domain by taking the difference

Second way: summing up all fluxes that leave/enter the domain, these are:

-advection of tracers across the OPBCs (I used the routine 
obcs_balance_flow.F as a template and added the multiplication by the 
pTracer value corresponding to first order upwind scheme, since I 'm 
setting OBCS_u1_adv_Tr=10 in data.obcs

-sedimentation of tracers, leaving the domain through the lower most wet 
grid cell (modification to DIC that I did myself)

-cumulating phosphorus content added to the domain due to resetting 
negative tracer values to zero

- RBCS contribution, by taking the difference in gPtracer in 
ptracers_apply_forcing.F after and before RBCS_ADD_TENDENCY is called 
(and then integrating it over volume of the corresponding grid cells and 
time step)

I manage to match the two different ways of calculating the total change 
in phosphorous in my domain only when I either close the domain with 
solid walls instead of open boundaries, or if I switch the advection and 
diffusion of ptracers of (and possibly other contribution that I'm not 
aware of, by either commenting out the call to Ptracers_integrate.f in 
thermodynamics.F or by setting gTracer to zero around line 380 in 
ptracers_integrate.F).


So I wonder whether I forgot some fluxes (is there diffusion across the 
open boundaries) or whether I calculate the fluxes across the open 
boundaries at a wrong time.  I tried to calculate them at the end of 
gchem_forcing_sep.F and at the end of ptracers_integrate.F just before 
OBCS_APPLY_PTRACER. (I assume that this call to OBCS_APPLY_PTRACER 
updates the ptracer with values of time-level n+1, since the advection 
across the OPBs for the current time level is done before).

I'm sorry for this very long and possibly confusing email. Any ideas on 
other fluxes contributing to ptracer budgets or hints on where in the 
routines to calculate the fluxes across the open boundaries will be 
greatly appreciated.

Kind regards,

Hannah

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20211009/0e2f4a44/attachment.html>