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

new CDEPS stream capability to read in forcing files needed by HAMOCC #405

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

Conversation

mvertens
Copy link
Contributor

@mvertens mvertens commented Sep 23, 2024

This PR has new stream functionality for forcing files.

  • ocn_stream_rivin.F90 - reads in riverin data without time coordinate
    • and maps it go blom grid using mapping file for custom mapping
      files for mosart->blom mapping
  • mod_blom_pio.F90
    • new pio capability needed by ocn_stream_rivin.f90
  • ocn_stream_dust.F90
    • uses inline CDEPS stream functionality to read in and interpolate
      iron dust deposition data
  • ocn_stream_ndep.F90
    • uses inline CDEPS sream functionality to read in and interpolatte
      nitrogen deposition data
    • in testing mode.
  • ocn_stream_chloro.F90
    • in testing mode
  • ocn_stream_swa.F90
    • in testing mode

@mvertens mvertens marked this pull request as draft September 23, 2024 08:48
@mvertens
Copy link
Contributor Author

@jmaerz @JorgSchwinger

  • Currently, CAM sends NDEP forcing to the mediator using prognostic output from WACCM OR reading it in from forcing files using CDEPS inline functionality. So BLOM should never try to read this data if it's running with CAM. I'll put this capability in this PR.

  • When BLOM runs with DATM it can use whatever ndep forcing it needs (or listen to DATM). Howver, the current ndep forcing files don't have the required time coordinates to be read in by CDEPS.

What is required is a time coordinate like the following:

	double time(time) ;
		time:long_name = "time" ;
		time:units = "days since 0001-01-01 00:00:00" ;
		time:calendar = "noleap" ;
		time:bounds = "time_bnds" ;
		time:cell_methods = "time: mean" ;

The follow files do not have this.

  <value blom_n_deposition="yes" blom_ndep_scenario="1850" ocn_grid="tnx2v1">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_1850_CMIP6_tnx2v1_20180321.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="2000" ocn_grid="tnx2v1">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_2000_CMIP6_tnx2v1_20200826.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="hist" ocn_grid="tnx2v1">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_185001-201412_tnx2v1_20190702.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="1850" ocn_grid="tnx1v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_1850_CMIP6_tnx1v4_20171106.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="2000" ocn_grid="tnx1v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_2000_CMIP6_tnx1v4_20200826.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="hist" ocn_grid="tnx1v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_185001-201412_tnx1v4_20180613.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="ssp"  ocn_grid="tnx1v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_201501-210012-${BLOM_NDEP_SCENARIO}_tnx1v4_20191112.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="1850" ocn_grid="tnx0.25v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_1850_CMIP6_tnx0.25v4_20190912.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="2000" ocn_grid="tnx0.25v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_2000_CMIP6_tnx0.25v4_20200826.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="hist" ocn_grid="tnx0.25v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_185001-201412_tnx0.25v4_20190705.nc</value>
  <value blom_n_deposition="yes" blom_ndep_scenario="1850" ocn_grid="tnx0.125v4">$DIN_LOC_ROOT/ocn/blom/bndcon/ndep_1850_CMIP6_tnx0.125v4_20221013.nc</value>

What I am currently doing is pointing to a forcing file that is used by CTSM:

$DIN_LOC_ROOT/lnd/clm2/ndepdata/fndep_clm_WACCM6_CMIP6piControl001_y21-50avg_1850monthly_0.95x1.25_c180802.nc

This file is used for all the current grids for pre-industrial runs using DATM.

@mvertens mvertens added the enhancement New feature or request label Sep 23, 2024
@mvertens mvertens added this to the NorESM2.5 - BLOM/iHAMOCC milestone Sep 23, 2024
stream_taxmode = 'cycle', &
stream_dtlimit = 1.0e30_r8, &
stream_tintalgo = 'linear', &
stream_name = 'Sea surface temperature ', &
Copy link
Collaborator

Choose a reason for hiding this comment

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

I syuspect this will be changed? - stream_name = Nitrogen deposition

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Right! Thanks. I'm just starting to test this today.

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 23, 2024

Dear @mvertens , many thanks for your efforts. The first thing that comes to my mind when initially browsing the changes: is there a good reason to use BLOM naming conventions for the grid size in iHAMOCC (i.e. exchanging iHAMOCC internally used kpie,kpje,.. by the BLOM idm and jdm)? From my side, I would remain with the iHAMOCC namings, but maybe @JorgSchwinger has a different opinion...

Wrt the input files for N-deposition: I feel it could be good to somehow organize the move from an input stream for NOy only to both NOy and NHx - can you tell me what your $DIN_LOC_ROOT is, so that I can have a look? As far as I understand, supporting currently still both, MCT and NUOPC means to provide input files also for different grids anyway - meaning that I will have to provide N-deposition data in a new manner, when moving to provide NOy and NHx input streams (and combining them in mo_apply_ndep). I try, but cannot promise to bring this further the next days due to time constrains (the scripts are already almost there, though - only the time variable might need to be added).

@JorgSchwinger
Copy link
Contributor

The follow files do not have this.

Regarding the time coordinates, my understanding was

  1. BLOM runs with NUOPC -> N-deposition comes from the NUOPC cap (either from CAM or DATM)
  2. BLOM runs with MCT -> we use the old reading of files inside HAMOCC

In both cases we don't need new/updated files. We might need to see if everything is available for DATM cases, but for now, if there is a pre-industrial N-deposition available that would be fine from my perspective.

@mvertens
Copy link
Contributor Author

@jmaerz - I am happy to change the names to kpie and kpje - but this code is no longer in the internal blom code and so all I will do is end up making a copy of these names in the nuopc cap.
i.e. kpie - idm, kpje = jdm.
If both of you think this will add clarity - I'll go ahead and make this change.

DIN_LOC_ROOT=/cluster/shared/noresm/inputdata. This already has NOy and NHx but with different names in the forcing data.

@mvertens
Copy link
Contributor Author

@JorgSchwinger - thanks for confirming what I was suggesting. Somehow I understood that you still wanted NDEP forcing read in even for nuopc. Using the datm functionality removes the need to have this in place. I just need to get this working in the cap. If @jmaerz is in agreement we won't need ndep stream functionality which using nuopc.

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 23, 2024

@JorgSchwinger , what we can do is to rewrite the mo_read_ndep just slightly - always providing/allocate both fields, nhx and noy, while only filling noy in case of the regular, old N-cycle - and then sum up noy and nhx in mo_apply_ndep for the old N-cycle. This would ease the whole procedure and no new files would need to be written, if I am, right now not mistaken (just to be checked, whether units, etc are correct in the other, newly interpolated input files that Mariana now sends through CDEPS).

@mvertens
Copy link
Contributor Author

@jmaerz - good point about units.
Looks like the BLOM current input is kmol N m-2 yr-1 (looking at ndep_1850_CMIP6_tnx1v4_20171106.nc)
The input from the dataset I am pointing to is kg/m2/s.
So this is confusing to me. If data is in a yr-1 units - what does it mean to do a monthly interpolation?

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 23, 2024

Well, kg/m2/s * (mol N/g) * (s/year) = kmol/m2/year would be the conversion. With mol-weight being 14g/mol, year=365*86400 (those brought my new processing closest to the 'old' N-deposition files). I am not sure, though, if we should make the move to i) year as defined by the model and ii) mol-weight as defined by CAM (14.00674 g/mol N)... - maybe for the latter even using the CAM chemical library (or use the value defined in mo_chemcon.F90 in hamocc).

To briefly summarize the discussion that I had during lunch with Jörg: We agreed on having by default both input-streams defined via nuopc, both, NOy and NHx - so similar to dust, etc. This makes it necessary to reformulate the N-deposition a bit in mo_apply_ndep.F90 (sum up NOy and NHx in case of the old N-cycle) and also the reading of input files in mo_read_ndep.F90 (always providing both NOy and NHx, while only reading NOy and NHx remains zero in case of the old files), which can both be easily done. This makes the checking of the nndep parameter obsolete - it always will be 2 (or your streams, equivalently).

@mvertens
Copy link
Contributor Author

mvertens commented Sep 23, 2024

Well, kg/m2/s * (mol N/g) * (s/year) = kmol/m2/year would be the conversion. With mol-weight being 14g/mol, year=365*86400 (those brought my new processing closest to the 'old' N-deposition files). I am not sure, though, if we should make the move to i) year as defined by the model and ii) mol-weight as defined by CAM (14.00674 g/mol N)... - maybe for the latter even using the CAM chemical library (or use the value defined in mo_chemcon.F90 in hamocc).

I see that the conversion is done in mo_read_ndep.F90.

To briefly summarize the discussion that I had during lunch with Jörg: We agreed on having by default both input-streams defined via nuopc, both, NOy and NHx - so similar to dust, etc. This makes it necessary to reformulate the N-deposition a bit in mo_apply_ndep.F90 (sum up NOy and NHx in case of the old N-cycle) and also the reading of input files in mo_read_ndep.F90 (always providing both NOy and NHx, while only reading NOy and NHx remains zero in case of the old files), which can both be easily done. This makes the checking of the nndep parameter obsolete - it always will be 2 (or your streams, equivalently).

I'm not sure I understand what you are suggesting. With NUOPC you always receive ndep (nhx and noy separately) from the mediator. This is either from DATM or CAM. Why do we need to read it in separately in the nuopc cap? And when would you trigger this? With the complexity of supporting both MCT and NUOPC and then asking for stream input from NDEP despite that you always receive it - I think the complexity of supporting all of this will get very confusing.

@mvertens
Copy link
Contributor Author

mvertens commented Sep 23, 2024

Here is one example of the complexity:

  1. With the updated NUOPC cap that I am putting together - you will always get NDEP from the mediator.

So it does not matter what HAMOCC_EXTNCYCLE or HAMOCC_NDEPC are set to.

 <entry id="HAMOCC_EXTNCYCLE">
    <type>logical</type>
    <valid_values>TRUE,FALSE</valid_values>
    <default_value>FALSE</default_value>
    <values>
       <value compset="_BGC.*%N2OC">TRUE</value>
       <value compset="_BGC.*%NH3C">TRUE</value>
       <value compset="_BGC.*%ATMNDEPC">TRUE</value>
    </values>
    <group>run_component_blom</group>
    <file>env_run.xml</file>
    <desc>Set namelist option to activate the extended nitrogen cycle code. Requires module ecosys</desc>
  </entry>

  <entry id="HAMOCC_ATMNDEPC">
    <type>logical</type>
    <valid_values>TRUE,FALSE</valid_values>
    <default_value>FALSE</default_value>
    <values>
       <value compset="_BGC.*%ATMNDEPC">TRUE</value>
    </values>
    <group>run_component_blom</group>
    <file>env_run.xml</file>
    <desc>Nitrogen deposition coupled from atmosphere. Requires module ecosys and extncycle</desc>
  </entry>
  1. In namelist_definition_blom.xml - if you sett the attribute modify_via_xml="" then you will ignore any value that is set in values. The value will come from HAMOCC_ATMNDEPC. So this does not make sense to me. The logic should be in setting HAMOCC_ATMNDEPC to TRUE. Am I missing something?
  <entry id="do_ndep_coupled" modify_via_xml="HAMOCC_ATMNDEPC">
    <type>logical</type>
    <category>bgcnml</category>
    <group>bgcnml</group>
    <values>
      <value>.false.</value>
      <value blom_n_deposition="yes" hamocc_atmndepc="yes" hamocc_extncycle="yes">.true.</value>
    </values>
    <desc>Switch to couple nitrogen deposition. Requires do_ndep.</desc>
  </entry>

  <entry id="do_n2onh3_coupled" modify_via_xml="HAMOCC_N2OC">
    <type>logical</type>
    <category>bgcnml</category>
    <group>bgcnml</group>
    <values>
      <value>.false.</value>
      <value hamocc_n2oc="yes" hamocc_extncycle="yes">.true.</value>
    </values>
    <desc>Switch to couple N2O and NH3 fluxes</desc>
  </entry>

I think that another call would be helpful again.

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 23, 2024

Dear @mvertens , I have feeling, we again misunderstand each other - sorry for this. All I want to say is that via nuopc, it would be good to receive both, NHx and NOy. Then, as you've done already for dust, if the nuopc stream is used (either from CAM or from DATM), ndep_stream feeds directly into mo_apply_ndep.F90. There however, it was previously assumed that the deposition holds both for the old N-cycle in one field (so NO3 deposition in iHAMOCC=NOy+NHx). Hence a summation is now needed there for the old N-cycle. In case of the new extended N-cycle, both fields are required.
In case that MCT is used, there is still the need to enter mo_read_ndep (like you implemented for dust). In the latter case, it was currently pre-defined, if one field or two fields are read in for N deposition - now we should always have a ndep_stream that has two components (NOy and NHx), while it either remains with filled NOy and zero NHx in case of the old Ncycle or both filled in the extended N cycle. The files in the bnd conditions for BLOM currently hold ndep=NOy+NHx (just for explanation - see also my intitial PR on simplifying things #398 - that gets obsolete with your new implementation, while the summation needs to be done in mo_apply_ndep).

@mvertens
Copy link
Contributor Author

@jmaerz - thanks for the clarification. That's very helpful. I totally misunderstood what you were saying.

call iniswa
if (.not. use_stream_swa) then
call iniswa
end if
Copy link
Contributor

Choose a reason for hiding this comment

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

I fear there has been a misunderstanding about the HAMOCC swa field? For the bromoform-option there has been a swa-file, which contains a pre-industrial swa-climatology (because bromoform production is parameterized using deviations from the swa rather than absolute values). This iniswa deals with how shortwave absorption is handled in BLOM and need to be called always.

The "normal" short-wave radiation field as far as I am aware always comes from CAM or DATM and is never read through a stream by BLOM.

The reading of the swa-climatology is not done consistently with the other inputs (the bromoform code has not been used a lot). The climatology file is read in in homocc_init and the just used inside hamocc.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@JorgSchwinger - thanks for clarifying. I should have seen this.

@@ -43,7 +43,7 @@ module mo_read_fedep

contains

subroutine ini_read_fedep(kpie,kpje,omask)
subroutine ini_read_fedep(omask)
Copy link
Contributor

Choose a reason for hiding this comment

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

I tend to agree with @jmaerz ... All HAMOCC routines are taking the indices as arguments, so why change this here?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@JorgSchwinger - I'll put the original back in. I was trying to make this consistent with the way streams was being handled.

@mvertens
Copy link
Contributor Author

mvertens commented Sep 23, 2024

As a sanity check (until I get the new dataset from @JorgSchwinger ) on the new river nutrient implementation I have used the file /cluster/shared/noresm/inputdata//ocn/blom/bndcon/river_nutrients_GNEWS2000_tnx1v4_20170820.nc
as the input and expect to see identical results using the compset NOINYOC and the resolution tnx1v4.
The good news is that I do. The left hand screen is data from river_nutrients_GNEWS2000_tnx1v4_20170820.nc and the right hand screen is what the ocn_stream_rivin.F90 produces.

Screenshot 2024-09-23 at 4 32 45 PM

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 23, 2024

@mvertens , thanks for testing! - yes, looks good to me as well :-)

@JorgSchwinger
Copy link
Contributor

Dear @mvertens, I have produced river nutrient and carbon inputs to be used by HAMOCC on the 0.5x0.5 runoff grid. I am not able to exactly reproduce the "original" files that are currently used for NorESM2. This is (most likely) because the scripts read data from xls-files and they are in matlab, which I don't have so I use octave... In general the deviations are about as expected under such circumstances. I also produced a file on the tnx1v4 grid for comparison.

The files are here (@TomasTorsvik can you make sure Mariana has access to NS2980K, takk?):

/projects/NS2980K/schwinger/NorESM_inputdata/RIVNUT/GNEWSregrid/river_nutrients_GNEWS2000c70_r05_20240921.nc
/projects/NS2980K/schwinger/NorESM_inputdata/RIVNUT/GNEWSregrid/river_nutrients_GNEWS2000c70_tnx1v4_20240921.nc

I also made a dust-file with correct time-coordinate, which you can find here:

/projects/NS2980K/schwinger/NorESM_inputdata/dust/NorESM2/dustdep_mhw2006_tnx1v4_20171107_newtime.nc

@TomasTorsvik
Copy link
Contributor

@JorgSchwinger , @mvertens , @jmaerz - I have given access to Mariana for the NS2980K storage project.

@mvertens
Copy link
Contributor Author

@JorgSchwinger - thanks so much for the files!

@mvertens
Copy link
Contributor Author

Using the following namelist in ocn_in and setting USE_NUOPC_RIVIN = .true.

&stream_rivin
  STREAM_RIVIN_DATA_FILENAME = '/cluster/shared/noresm/inputdata/ocn/blom/bndcon/river_nutrients_GNEWS2000c70_r05_20240921.nc'
  STREAM_RIVIN_MAP_FILENAME = '/cluster/shared/noresm/inputdata/cpl/cpl6/map_r05_to_tnx1v4_e1000r300_170609.nc'
  STREAM_RIVIN_MESH_FILENAME = '/cluster/shared/noresm/inputdata/share/meshes/r05_nomask_c110308_ESMFmesh.nc'
/

I get the following for DIC:

Screenshot 2024-09-25 at 10 34 29 AM


dust(:,:) = dustflx(:,:,kplmon)

if (debug) then
Copy link
Collaborator

Choose a reason for hiding this comment

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

Hi @mvertens , do you want to keep the debug flag here at the end or is it something temporary? If the first, I wonder, if we should use our common hamocc debug flag(s) instead.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's up to you. Do you always want to see the netcdf file when the model starts?

@@ -172,32 +172,28 @@ subroutine ini_read_ndep(kpie,kpje)
end subroutine ini_read_ndep


subroutine get_ndep(kpie,kpje,kbnd,kplyear,kplmon,omask,ndep,patmnhxdep,patmnoydep)
subroutine get_ndep(kplyear, kplmon, omask, ndep, patmnhxdep, patmnoydep)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Here the same wrt kpie,...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done. Thanks for catching this.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I wasn't sure, if you had already reverted back or not - there were/are a few occasions, were it was introduced.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done.


!***********************************************************************************************
! Read and return CMIP6 n-deposition data for a given month or use atmosphere input
!
! S. Gao *Gfi, Bergen* 19.08.2017
!***********************************************************************************************

use mod_xc, only: mnproc
use mod_xc, only: idm, jdm, nbdy, mnproc
Copy link
Collaborator

Choose a reason for hiding this comment

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

Again about kpie, ...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done.

@@ -262,30 +259,29 @@ subroutine ini_read_oafx(kpie,kpje,pdlxp,pdlyp,pglat,omask)
end subroutine ini_read_oafx


subroutine get_oafx(kpie,kpje,kplyear,kplmon,omask,oafx)
subroutine get_oafx(kplyear, kplmon, omask, oafx)
Copy link
Collaborator

Choose a reason for hiding this comment

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

kpie,...

@mvertens
Copy link
Contributor Author

@jmaerz - I put the debug in order to trigger generating the netcdf file showing the forcing data. I'm happy to remove it - then you'll always get the file. Is that what you would want?

@jmaerz
Copy link
Collaborator

jmaerz commented Sep 25, 2024

I suspect that we at the end would want to have some general output of the input fields (since it will be now regridded internally) - for now, I was just wondering. I would expect that we do the output then at a different place. All fine.

@@ -172,15 +172,15 @@ subroutine ini_read_ndep(kpie,kpje)
end subroutine ini_read_ndep


subroutine get_ndep(kpie,kpje,kbnd,kplyear,kplmon,omask,ndep,patmnhxdep,patmnoydep)
subroutine get_ndep(idm, jdm, kplyear, kplmon, omask, ndep, patmnhxdep, patmnoydep)
Copy link
Collaborator

Choose a reason for hiding this comment

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

kpie,...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
Status: Todo
Development

Successfully merging this pull request may close these issues.

4 participants