Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Read BGC input files through pkg/gchem #876

Draft
wants to merge 20 commits into
base: master
Choose a base branch
from
Draft

Conversation

averdy
Copy link
Collaborator

@averdy averdy commented Oct 15, 2024

What changes does this PR introduce?

Feature: BGC input files are handled through pkg/gchem. A more general approach than what was initially proposed for pkg/bling ( PR #841 ), as it can be used for other BGC models. Addresses issue #138

What is the current behaviour?

BGC forcing (surface iron flux, atmospheric pCO2, wind speed, atmospheric pressure, ice) and fields that are not prognostic variables (surface silica) are specified in namelist files data.dic, data.bling, or data.darwin and read in the package-specific *_fields_load.F routine.

What is the new behaviour

Those fields are now read in gchem_fields_load.F and handled through data.gchem. It is done using pkg/exf routines if useEXF (or maybe useGCHEM_EXF?) is true; otherwise it is done using either the current method (READ_REC_XY_RS).

Does this PR introduce a breaking change?

For pkg/bling, the only change required is switching BGC forcing fields info to data.gchem. For other packages, no breaking change - by default, forcing fields read with gchem will not be passed on to the package.

Other information:

  • For winds and atmos pressure, need to check compatibility of exf fields (wspeed and apressure) with BGC code
  • Currently doesn't run with #define USE_EXF_INTERPOLATION
  • Eventually we want to read in 3D silica fields

Suggested addition to tag-index

@mjlosch
Copy link
Member

mjlosch commented Oct 15, 2024

Hi @averdy
thanks you for making a great start here. Could you please modify this PR, so that it is easier to compare, e.g.

  • the directory gchem should contain all the fields necessary to compile the model (I guess it does), could you please remove gchem/gchem_orig (we don't need the originals here, git takes care of that).
  • the new verification experiment input.exf contains duplicates of input and input_ad. I assume that it is a fwd experiment, so the ad-specific file are maybe not needed at all, and the rest of the duplicates should be (automatically) linked from input, so they can be removed, too, e.g. data.gmredi. It makes it easier to concentrate on the relevant changes.
  • If you add an output.exf.txt to results, it would be easier to test this

I'll have a close look, once this is done, but for naming I suggest that the field names need not necessarily include "exf", because one can image reading the same fields via the "old" method. For a cleaner separation, we could think of having a separate routine gchem_exf_fields_load or gchem_exf_getffields.F that is called from (or before) gchem_fields_load.

@jm-c
Copy link
Member

jm-c commented Oct 15, 2024

@averdy Thanks for starting this. I suggest to wait a little bit before any new updates here to give us some time to provide feedback.

@mjlosch
Copy link
Member

mjlosch commented Oct 15, 2024

I needed to have a data with

#tauThetaClimRelax = 5184000.,
#tauSaltClimRelax  = 7776000.,

commented out, but with these changes:

diff --git a/pkg/gchem/gchem_init_vari.F b/pkg/gchem/gchem_init_vari.F
index a877b1c53..d5926e406 100644
--- a/pkg/gchem/gchem_init_vari.F
+++ b/pkg/gchem/gchem_init_vari.F
@@ -169,8 +169,8 @@ C-- Initialize Geo-Chemistry forcing fields
         _EXCH_XY_RL( gchemSi, myThid )
         _EXCH_XY_RL( gchemPAR, myThid )
         _EXCH_XY_RL( gchemFe, myThid )
-        _EXCH_XY_RL( gchemIce, myThid )
-        _EXCH_XY_RL( gchemWind, myThid )
+c       _EXCH_XY_RL( gchemIce, myThid )
+c       _EXCH_XY_RL( gchemWind, myThid )
         _EXCH_XY_RL( gchemapCO2, myThid )
 
 #endif /* ALLOW_EXF */
diff --git a/pkg/gchem/gchem_readparms.F b/pkg/gchem/gchem_readparms.F
index 8968149f4..bff3ce516 100644
--- a/pkg/gchem/gchem_readparms.F
+++ b/pkg/gchem/gchem_readparms.F
@@ -102,7 +102,7 @@ C- Open and read the data.gchem file
      I                   myThid )
       READ(UNIT=iUnit,NML=GCHEM_PARM01)
 #ifdef ALLOW_EXF
-      CALL GCHEM_EXF_READPARMS( iUnit, myThid )
+      IF (useEXF) CALL GCHEM_EXF_READPARMS( iUnit, myThid )
 #endif /* ALLOW_EXF */
       WRITE(msgBuf,'(A)')
      &  ' GCHEM_READPARMS: finished reading data.gchem'
diff --git a/verification/global_oce_biogeo_bling/input.exf/data.pkg b/verification/global_oce_biogeo_bling/input.exf/data.pkg
index 2564a8df7..04701bb30 100644
--- a/verification/global_oce_biogeo_bling/input.exf/data.pkg
+++ b/verification/global_oce_biogeo_bling/input.exf/data.pkg
@@ -3,7 +3,7 @@
  useCAL=.TRUE.,
  useEXF=.TRUE.,
  useGMRedi=.TRUE.,
- usePROFILES=.TRUE.,
+#usePROFILES=.TRUE.,
  usePTRACERS=.TRUE.,
  useGCHEM=.TRUE.,
  useDiagnostics=.TRUE.,

I can now run input.exf, but the default fwd experiment fails.
@averdy Did you try that?

@averdy
Copy link
Collaborator Author

averdy commented Oct 15, 2024

@mjlosch I'm having computing issues and can't run the model right now, but I'll look into it as soon as possible! I'm sharing this as a draft mostly to get feedback on the general structure. I'm not attached to any of it at this point!

averdy and others added 3 commits October 15, 2024 09:46
avoids using unitialised filenames that make reference experiment
global_oce_biogeo_bling fail (and it is better style)
Copy link
Member

@mjlosch mjlosch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that this is great. I have a few inline comments and suggestions.

I would even go a step further and make pkg/bling (and pkg/dic in a similar way) entirely agnostic of where the forcing fields come from. A gchem_fields_load loads all all forcing data provided by data.gchem into fields called gchem_wspeed, etc, and it is here, where we decide if we want to use the ext-method or the original simple method. bling_fields_load would then just copy to forcing fields local to bling (or dic), or uses the gchem-fields explicitly (i.e. bling_fields_load would no longer be necessary). Each BGC package called from gchem can use the fields or decide not to use them and have their own method of loading them.

pkg/bling/bling_fields_load.F Outdated Show resolved Hide resolved
pkg/bling/bling_fields_load.F Outdated Show resolved Hide resolved
STOP 'ABNORMAL END: S/R INSOL: 2-D output not implemented'
ENDIF

C Find day (****NOTE for year starting in winter*****)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This goes a little further, but this bit is basically a copy of dic/insol.F, isn't it? Would it make sense to have a routine pkg/gchem/gchem_insol.F (or a similar name) that then handles the computation of irr_surf if irr_surf is not read from file?
Or even have this routine pkg/gchem/gchem_insol.F return some irr_surf that is either computed from these astronomical parameters or retrieved from a file, either by exf or not, it could even include the QSW_underice part.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mjlosch Yes, it's a copy of dic/insol.F and I agree that it would make sense to have a common routine

@mjlosch
Copy link
Member

mjlosch commented Oct 16, 2024

I would even go a step further and make pkg/bling (and pkg/dic in a similar way) entirely agnostic of where the forcing fields come from.

I have a version in my local repo, where I have done this. It gives exactly the same solution (via testreport) as the current code.
To show this to you without modifying this branch I am posting a patch that you can apply to your repo by

git apply gchem-bling-exf.patch

gchem-bling-exf.patch

It's not complete (e.g. some of the runtime parameters such as BLING_pCO2 need to be transferred to gchem, etc. but the main structure is what I meant in my earlier description. If we were to go in this direction, one would need to sort out the namelists (no clear separation between filenames with and without exf), etc.

averdy and others added 4 commits October 22, 2024 09:13
- bling_fields_load.F only contains a copy and possibly a unit
  conversion
- replace exf_ prefix by gchem_ (not fully consistent yet)
- always call gchem_exf_readparams
- gchem_fields_load.F contains all reading of forcing files, with and
  without exf; logic may not yet be optimal when to read with pkg/exf
- some sensible initialisation of fields in gchem_init_fixed.F is
  necessary to main results of reference experiment
- adjust namelists of reference experiment(s)
- add preliminary reference output for testing
- this does not yet work with TAF AD because variable names have
  changed and need adjustment; postponed to when naming changes are complete
tutorial_cfc_example and tutorial_global_oce_biogeo now compile and run
@mjlosch
Copy link
Member

mjlosch commented Oct 29, 2024

I did not pay to much attention to other verification experiments than global_oce_biogeo_bling, so that there were some problem, when bling was not used. The fixes in 82ea1c5 addressed this, but we may want to think about a different organisation of gchem_readparms.F and gchem_exf_readparms.F.

I note that none of this works with TAF, because variable names have changed and variable have moved. We will sort this out later, when we have converged on a complete set of names of the new gchem forcing variables.

A few thoughts:

  • We could use gchem_ or gc_ to prepend variables where necessary. I prefer gchem_ (pkg-name).
  • Many of the new variables will only be used in gchem_fields_load.F, especially gSI0/1 etc., so we may want keep them in a separate header file from the actual forcing field gchemSi (or maybe gchem_Si).
  • The name GCHEM_EXF.h may also not be a good choice, because these fields are also used if pkg/exf is not used.
  • I suggest to move the forcing fields gchemSi (and all others) to GCHEM.h or a new file GCHEM_FORCING.h and rename GCHEM_EXF.h to GCHEM_FIELDS.F or similar for the remaining gSi0/1 etc.

mjlosch and others added 4 commits October 30, 2024 12:13
- rename forcing fields to have prefix "gchem" (for now without _)
- forcing file names in general gchem namelist, maybe it is even
  better to move all of gchem_exf_readparms.F to gchem_readparms.F,
  including the forcing namelist?
- foward testreport passes (at least on my computer)
@mjlosch
Copy link
Member

mjlosch commented Nov 7, 2024

@averdy the latest commit breaks some of the verification experiments with gchem that do now use pkg/exf (and don't compile it either). Maybe some #ifdef ALLOW_EXF are missing?

@averdy
Copy link
Collaborator Author

averdy commented Nov 7, 2024

@mjlosch Fixed, thanks!

@mjlosch
Copy link
Member

mjlosch commented Nov 8, 2024

@averdy not sure why, but 2495e36 and a63c177 change the results of verification/global_oce_biogeo_bling (and they don't change back afterwards). Do you have any idea why? It would be good to fix this.

@mjlosch
Copy link
Member

mjlosch commented Nov 8, 2024

The largest changes are from a63c177 where initial conditions are changed. I think we should stick to the old initial conditions to be able to reproduce the old results for now, so that we don't get confused. In a separate PR we can still change these defaults if really necessary.
The remaining changes are from 2495e36. They are not big, but still it would be good to understand why this happens (and potentially fix it):

default 10  ----T-----  ----S-----  ----U-----  ----V-----  --PTR 01--  --PTR 02--  --PTR 03--  --PTR 04--  --PTR 05--
G D M    c        m  s        m  s        m  s        m  s        m  s        m  s        m  s        m  s        m  s
e p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
n n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16  8 16 10> 7<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 FAIL  global_oce_biogeo_bling
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16  8 16 10> 7<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 FAIL  global_oce_biogeo_bling.exf

@mjlosch
Copy link
Member

mjlosch commented Nov 11, 2024

The remaining changes are from 2495e36. They are not big, but still it would be good to understand why this happens (and potentially fix it):

It has to do with the initialisation of Silica, but I don't understand the context ...

@jahn
Copy link
Member

jahn commented Nov 12, 2024

Not sure I understand why gchem_pco2 is needed. The exf package already has a mechanism for setting a constant value. Why not use apCO2const? If I understand the proposed changes, this is a runtime parameter but ignored. This is kind of unexpected. Same for the hard-coded default values for the other fields set in gchem_init_varia.F. Apologies if I am missing something, just trying to catch up.

@mjlosch
Copy link
Member

mjlosch commented Nov 12, 2024

There was bling_pCO2, so I thought I'd replace it. This is for the case that pkg/exf is not defined.

I totally see that this is redundant, and we can use apCO2const instead (but then need to reorganise name lists and common blocks). For this we'd probably need to merge gchem_readparms.F and gchem_exf_readparms.F, or at least have the GCHEM_EXF.h parameters defined all the time (and included in gchem_readparms.F).

The hard coded default values are only here to be able to reproduce old results (which already does not work anymore). This is now modelled a bit after the bling and dic-verification experiments, because that's what we have. There is no problem in modifying the defaults as soon as there's a consensus what they should be (probably not zero in some cases).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants