[MITgcm-support] Simulation with varying atmospheric pCO2

Wed Aug 18 12:21:02 EDT 2021

Dear Jonathan,

Thank you for your reply. I’m not sure to understand, what do you mean 
by log call ? I don't see conditions to limit a certain value in the 
files you pointed out.

I started looking on your suggestions and what I can say for now is:

(1) I copy the two files for ocean at the end of this message.

(2a) By only checking available outputs every half-an-hour, it seems the 
atmospheric pCO2 concentration becomes equals to zero first. 
Unfortunately, as this is a homogeneous value for the whole atmosphere I 
cannot identify a local point where the problem occur. But following 
your suggestion (4) it seems the problem do not directly coming from here.

(2b) For the moment, I haven’t observe the problem but both atmospheric 
and oceanic pCO2 concentration are increasing. My simulation is still 
running so I will see if this solve the problem after spin-up time or not.

(3) Yes, I tried and it also turn to the same problem.

(4) Thank you for this suggestion. Looking on the additional prints, 
especially those related to dic_atmos.F file, leads me to believe that 
the problem first occur with the fluxCO2 variable in ocean - which is 
shared with atmosphere. Here I copy the relevant lines of debug output.

QQ atmos C, total, pCo2 NaN NaN
QQ total C, current, old, diff NaN 1.3077058455813389E+018 NaN
QQ ocean C, current, old, diff NaN 1.2264261974040456E+018 NaN
QQ air-sea flux, addition diff NaN NaN
DEBUG_STATS_RL: atpco2 (DIC_ATMOS) min= 0.00000000000000E+00
DEBUG_STATS_RL: atpco2 (DIC_ATMOS) max= 0.00000000000000E+00
DEBUG_STATS_RL: atpco2 (DIC_ATMOS) mean= 0.00000000000000E+00
DEBUG_STATS_RL: atpco2 (DIC_ATMOS) S.D.= 0.00000000000000E+00
DEBUG_STATS_RL: co2atmos (DIC_ATMOS) min= 0.00000000000000E+00
DEBUG_STATS_RL: co2atmos (DIC_ATMOS) max= 0.00000000000000E+00
DEBUG_STATS_RL: co2atmos (DIC_ATMOS) mean= 0.00000000000000E+00
DEBUG_STATS_RL: co2atmos (DIC_ATMOS) S.D.= 0.00000000000000E+00
DEBUG_WRITE: atm_pCO2 val= NaN
DEBUG_WRITE: dic_pCO2 val= 3.20000000000000E-04
DEBUG_WRITE: tile_flux val= 1.07455380469069E+09
DEBUG_WRITE: total_atmos_carbon val= NaN
DEBUG_WRITE: total_flux val= NaN
DEBUG_STATS_RL: dTtracerLev min= -2.00000000000000E+03
DEBUG_STATS_RL: dTtracerLev max= 4.65311520000000E+09
DEBUG_STATS_RL: dTtracerLev mean= 2.39551968956715E+07
DEBUG_STATS_RL: dTtracerLev S.D.= 3.18669211791863E+08

There is points with fluxCO2 around -6.10⁻⁴ (so two order of magnitude 
smaller than elsewhere) at the limit between two faces (note there is no 
ice into this simulation). If I fix a lower limit for fluxCO2 then there 
is not anymore appearance of zeros, and atmospheric pCO2 stay in 
reasonable range. But, the oceanic pCO2 is largely increasing (up to 
8000 ppm!) and there is still “NaN” in standard output related to 
oceanic carbon. I put below the corresponding lines of standard outputs:

QQ atmos C, total, pCo2 81393757432549728. 4.5985173690706060E-004
QQ total C, current, old, diff NaN 1.3077629004370186E+018 NaN
QQ ocean C, current, old, diff NaN 1.2264261974040456E+018 NaN
QQ air-sea flux, addition diff -57054399576804.773 NaN

Is this flux can really be the initial source of the problem or am I 
wrong? Is make sense to limit it? I would have a look on it computation 
to better understand the problem and solve it but I can't found it, do 
you not in which file I should look?

Best regards,

Charline

Here are the files:

data.dic:

# DIC parameters
&ABIOTIC_PARMS
selectPHsolver=1
&

&BIOTIC_PARMS
alphaUniform = 0.97e-10,
rainRatioUniform = 5.e-2,
KRemin = 0.95,
&

&DIC_FORCING
dic_int1 = 3,
dic_pCO2 = 320.E-6,
DIC_iceFile=' ',
DIC_windFile=' ',
DIC_silicaFile=' ',
DIC_atmospFile=' ',
&

data.cpl:

# Coupling package parameters, OCN component:
# cpl_earlyExpImpCall :: call coupler early in the time stepping call 
sequence
# useImportHFlx :: True => use the Imported HeatFlux from couler
# useImportFW :: True => use the Imported Fresh Water flux fr cpl
# useImportTau :: True => use the Imported Wind-Stress from couler
# useImportSLP :: True => use the Imported Sea-level Atmos. Pressure
# useImportSIce :: True => use the Imported Sea-Ice loading
# cpl_taveFreq :: Frequency^-1 for time-Aver. output (s)
&CPL_OCN_PARAM
# cpl_earlyExpImpCall=.FALSE.,
# useImportHFlx=.FALSE.,
# useImportFW =.FALSE.,
useImportTau =.TRUE.,
useImportSLP =.FALSE.,
# useImportSIce=.FALSE.,
# 1 month average
# cpl_taveFreq=2592000.,
# 10 year average
cpl_taveFreq=31104000000.,
&

--
Charline RAGON
PhD Student
Group of Applied Physics & Institute for Environmental Sciences
University of Geneva, Switzerland

> Dear Charline,
>
> Thanks for your email, I?ll reply here for posterity.
>
> I?m not as familiar with the coupled set up as I would like, however we should be able to diagnose what?s going on. In my experience, abrupt NaNs in output like this are usually the result of illegal values in a log call. There?s one in ?pkg/aim_v23/phy_driver.F? and another in the pH solver for oceanic pCO2 in ?pkg/dic/carbon_chem.F? that I would target. I can see how the former could be the issue if aim_pco2 is invalid (negative or zero), which could happen if there is an issue in ?pkg/aim_v23/aim_do_co2.F?...I think you?d get warnings/errors or immediate crashes if your variable choices were an issue, and also the pCO2 appears to be increasing in your run with ?aim_select_pcO2 = 3?, not trending to zero,  but it wouldn?t hurt to check that the DIC package is indeed talking to the AIM package (for example, do you have "#define COMPONENT_MODULE? in your CPP_EEOPTIONS.h file?) On the other hand, the latter log-call is often at the root of many issues like this so I am wondering if th
>   e problem instead lies in the DIC package itself.
>
> So, a couple of things to test:
> (1) Can you post your "data.cpl" and "data.dic" files?
>
> (2a) Have a look at the ocean DIC/surface pCO2 and see if there are NaNs propogating there too. If you have fine enough output frequency around the crash you observe, you may be able to narrow down where geographically the issue is coming from (possibly a sea-ice region?).
>
> (2b) If you are using the default pCO2 solver (which can be bit finicky), I would recompile with "#define CARBONCHEM_SOLVESAPHE" in ?DIC_OPTIONS.h? and set "selectPHsolver=1? in the ?ABIOTIC_PARMS? section of ?data.dic?. This uses a more robust (although possibly more expensive) pH/pCO2 solver.
>
> (3) Have you tried restarting the model with variable absorption coefficient from the end of the run with constant absorption coefficient? I think sometimes there are too many moving parts which causes hard to pin-down erratic behavior as the model adjusts to the initial conditions.
>
> (4) If the above doesnt work you could use the debug framework (#define ALLOW_DEBUG in CPP_OPTIONS.h, set debugMode = .TRUE. in ?eedata? file, and set "debugLevel=5? in PARM01 section of your ?data? file ? lots of messages), plus adding a few print statements for identifying where the NaN first appears, such as maybe print out the relevant variables in ?pkg/aim_v23/aim_do_co2.F? and ?pkg/aim_v23/phy_driver.F?, as well as in ?pkg/dic/dic_atmos.F? and ?hnew?, ?gam? and ?stuff? variables in ?pkg/dic/carbon_chem.F? subroutine "CALC_PCO2_APPROX?.
>
> Hope that helps!
>
> Jonathan
> ___________________________________________________________________________
> Dr. Jonathan M. Lauderdale
> Research Scientist
> Department of Earth, Atmosphere and Planetary Sciences
> Massachusetts Institute of Technology
> 77 Massachusetts Avenue (54-1518)
> Cambridge, MA 02139, USA
> Office: +1 617 324 1568
> Cell   : +1 617 304 5661
> Email:mailto:jml1 at mit.edu  <mailto:jml1 at mit.edu>
> Web:http://paocweb.mit.edu/people/jml1  <http://paocweb.mit.edu/people/jml1>
> Twitter:https://twitter.com/jon_lauderdale
>   <https://twitter.com/jon_lauderdale>Blog:https://seamanticscience.wordpress.com/  <https://seamanticscience.wordpress.com/>
> Git:https://github.com/seamanticscience  <https://github.com/seamanticscience>
> ___________________________________________________________________________
>
>> On Jul 19, 2021, at 8:09 AM, Charline Ragon<Charline.Ragon at unige.ch>  wrote:
>>
>> Dear all,
>>
>> I'm performing simulations in a coupled atmosphere, ocean and sea-ice configuration, using the biogeochemical packages GCHEM, DIC and PTRACERS.
>>
>> I would like to allow a CO2-flux between atmosphere and ocean but I have troubles as all diagnostics are suddenly moving to zero (NaN are appearing in standard outputs). I set parameters dic_int1 = 3 and aim_select_pCO2 = 3, but if I change the value of aim_select_pCO2, this problem is not appearing (at least, I haven't observed it).
>>
>> I had a look on, and attached to this e-mail, the evolution of atmospheric pCO2 in three similar simulations where I only change the aim_select_pCO2 value to 1, 2 or 3. These runs start from nIter0=0 and tracers are initialized by files. Note that, in my case aim_pCO2 = aim_ref_pCO2 = 320 ppm so it seems the absorption in CO2-band is constant and equals 4 also when aim_select_pCO2 = 1 (so, similar to aim_select_pCO2 = 0).
>>
>> As I understand it, when the absorption in CO2-band is prescribed the simulation no exhibits the problem even if partial pressure of CO2 is largely increasing in both atmosphere and ocean. So, the zero-appearing should be related to this coefficient but I'm unable to understand exactly what the problem is nor to solve it.
>>
>> I'm wondering if someone know more about this, and maybe have suggestions for me ?
>>
>> Thanks in advance for your help!
>>
>> Best,
>> Charline
>>
>>
>> --
>> Charline RAGON
>> PhD Student
>> Group of Applied Physics & Institute for Environmental Sciences
>> University of Geneva, Switzerland

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210818/2c7e36a1/attachment-0001.html>