[MITgcm-support] Simulation with varying atmospheric pCO2
Charline Ragon
Charline.Ragon at unige.ch
Wed Sep 15 12:55:49 EDT 2021
Dear Jonathan,
First of all, apologize for time reply.
I printed out «hnew» values just before the problem occur, and there is
no negative value. I also checked «stuff» and «gamm» variables, but
values are not drastically changing from previous iteration so I guess
the problem not coming from here.
DEBUG_WRITE: HNEW val= 1.90372358090375E-08
DEBUG_WRITE: STUFF val= 2.93023895674188E-14
DEBUG_WRITE: GAMM val= 9.27725921370696E-01
The test with Munhoven solvesaphe routine is still running: I do not
observe the problem of NaN but both atmospheric and oceanic pCO2 are
continuously increasing, as you can see in the attached figure, which
looks strange for me. I was expected two opposite behavior (one
increasing, the other decreasing) until have an almost similar pressure
of CO2 into the two components. It seems there is a generation of CO2
but, if I run the same simulation (especially, same initial conditions)
with constant atmospheric pCO2, the oceanic reservoir is not increasing
as much, and stabilize with value not so far of atmospheric one. Do you
think is this a normal behavior ?
I also tested adding the silicate forcing file used in
tutorial_global_oce_biogeo, but interpolated on my grid. The problem is
also occurring, at the same time step.
Regards,
Charline
--
Charline RAGON
PhD Student
Group of Applied Physics & Institute for Environmental Sciences
University of Geneva, Switzerland
> Dear Charline,
>
>
> Thanks for those additional information and update.
>
>
> I agree, it looks to me as though there is a problem occurring in your CO2 flux calculation, but if you are imposing the CO2 flux and still getting NaNs in the ocean carbon system then I suspect that there is a different underlying problem.
>
>
> The CO2 flux evaluation occurs in dic_surfforcing.F. Here, the model needs to calculate surface ocean pCO2 (either SR calc_pco2_approx in carbon_chem.F, or calc_pco2_solvesaphe in dic_solvesaphe.F, but let's focus on calc_pco2_approx), which requires solving for ocean pH via evaluation of hydrogen ion concentration. It's uncommon, but possible, for the calculation of [H+] ions (lines 438-443 in carbon_chem.F) to produce negative values, which results in NaNs when pH is evaluated (-log10([H+])) on line 448.
>
>
> To confirm negative [H+], I would put a "DEBUG_WRITE" statement in carbon_chem.F to print out the variable "hnew", or at least warn/STOP at negative values. There may be lots of output, so you could pick up from a point just before the model crashes.
>
>
> I'm hoping that the Munhoven (2013) solvesaphe routine will remedy this because it is supposed to be more robust. However, I notice that in your data.dic, your "DIC_silicaFile" variable is empty. Silicate concentration has a minor, but important, contribution to ocean pH, so I would recommend supplying an input file for that, otherwise the model assumes a constant, spatially invariant value which could be at the root this unexpected behaviour.
>
>
> Usually, we have supplied a spatially-varying monthly (annually repeating) field of concentrations from climatology (but it should be consistent with your "externForcingPeriod" and "externForcingCycle", if you set these, or you can force a 12 month repeating cycle by setting "DIC_forcingPeriod" and "DIC_forcingCycle" instead). Try this with the default pCO2/pH solver and see if that helps?
>
>
> Best of luck,
>
>
> Jonathan
>
> ___________________________________________________________________________
> Dr. Jonathan M. Lauderdale
> Research Scientist
> Department of Earth, Atmosphere and Planetary Sciences
> Massachusetts Institute of Technology
> 77 Massachusetts Avenue
> Cambridge, MA 02139, USA
> Office: +1 617 324 3401
> Cell : +1 617 304 5661
> Email: mailto:jml1 at mit.edu
> Web: http://paocweb.mit.edu/people/jml1
> Twitter: https://twitter.com/jon_lauderdale
> Git: https://github.com/seamanticscience
> ___________________________________________________________________________
> ________________________________
> From: MITgcm-support <mitgcm-support-bounces at mitgcm.org> on behalf of Charline Ragon <Charline.Ragon at unige.ch>
> Sent: 18 August 2021 12:21
> To: mitgcm-support at mitgcm.org
> Subject: Re: [MITgcm-support] Simulation with varying atmospheric pCO2
>
>
> 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><mailto:jml1 at mit.edu>
> Web: http://paocweb.mit.edu/people/jml1 <http://paocweb.mit.edu/people/jml1><http://paocweb.mit.edu/people/jml1>
> Twitter: https://twitter.com/jon_lauderdale
> <https://twitter.com/jon_lauderdale><https://twitter.com/jon_lauderdale>Blog: https://seamanticscience.wordpress.com/ <https://seamanticscience.wordpress.com/><https://seamanticscience.wordpress.com/>
> Git: https://github.com/seamanticscience <https://github.com/seamanticscience><https://github.com/seamanticscience>
> ___________________________________________________________________________
>
>
>
> On Jul 19, 2021, at 8:09 AM, Charline Ragon <Charline.Ragon at unige.ch><mailto: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 --------------
A non-text attachment was scrubbed...
Name: pCO2.png
Type: image/png
Size: 19928 bytes
Desc: not available
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210915/9868979a/attachment-0001.png>
More information about the MITgcm-support
mailing list