<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Dear MITgcm users and developers, <br>
</p>
<p>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.</p>
<p>1.) Time stepping of the DIC-biogeochemistry model:</p>
<p>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. </p>
<p>In the documentation I found the following sentence regarding
multi-dimensional advection (section 2.17.2.4)</p>
<p>"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 nth<span
class="math notranslate nohighlight"></span> time-level and not
the updated n+3/3<span class="math notranslate nohighlight"></span></p>
<p> quantities".</p>
<p><br>
</p>
<p>The full tracer equation should be sth like <br>
</p>
<p>Tr^(n+1)=Tr^n+GTr_adv^(n+1/2)*dt+GTr^n_diff*dt +GTr_bio^n *dt</p>
<p>the biological sink and source tendencies (GTr_bio) should then
be calculated based on TR^n<br>
</p>
<p>However, I believe the Tracer when entering subroutine
dic_biotic_forcing.F is</p>
<p>Tr*:= Tr^n+GTr_adv^(n+1/2)*dt+GTr_diff^n*dt</p>
<p>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?</p>
<p>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</p>
<p><br>
</p>
<p>2.) OBCS and ptracers</p>
<p> I calculated the total change of phosphorous in my domain in two
different ways:<br>
</p>
<p>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 <br>
</p>
<p>Second way: summing up all fluxes that leave/enter the domain,
these are:</p>
<p>-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 </p>
<p>-sedimentation of tracers, leaving the domain through the lower
most wet grid cell (modification to DIC that I did myself)</p>
<p>-cumulating phosphorus content added to the domain due to
resetting negative tracer values to zero</p>
<p>- 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)</p>
<p>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).<br>
</p>
<p><br>
</p>
<p>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).<br>
</p>
<p>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.</p>
<p>Kind regards,</p>
<p>Hannah</p>
</body>
</html>