Version 13 (modified by frrh, 4 years ago) (diff)

Catering for MEDUSA coupling via OASIS3-MCT


We need to cater for MEDUSA coupling at NEMO vn3.6 via the existing NEMO OASIS3-MCT interface and the Met Office UM code.

This work is initially primarily for the purposes of the UKESM project, however ultimately it is considered that the MEDUSA coupling interface would become a generally available option to any configuration which chooses to employ it.


The MEDUSA code itself is held in a development branch, created and maintained at NOC named, at time of writing, dev_r5518_NOC_MEDUSA_Stable This code potentially provides output fields to a coupler and receives input fields from a coupler to allow it to couple with an atmosphere model featuring dust and vegetation schemes.

Currently the specific coupling fields are expected to be:

  • Output:
    • CO2 Flux (MEDUSA field f_co2flux2d)
    • DMS (MEDUSA field dms_surf2d)
  • Input:
    • CO2 partial pressure (MEDUSA field f_pco2a. This is expected to come from the atmosphere and in the case of the UM specifically from the "TRIFFID" scheme)
    • Dust (MEDUSA field dust. This is expected to come from the atmosphere and in the case of the UM specifically from some special code which is as yet undeveloped!)

The code development needs to be compatible with the MO GO6 ocean model configuration. This configuration currently features approximately 11 separate code branches which makes further development of additional code potentially problematic in terms of a high risk of code conflicts.

It also needs to be compatible with the MO GC3 coupled model configuration which adds a further 4 NEMO branches.

While the bulk of the MEDUSA code is self contained under the TOP_SRC directory and thus less prone to causing conflicts, the coupling interface code is a very high risk in terms of conflicts since it needs to modify sbccpl.F90 and perhaps other high level coupling routines.

Oddly we find that the conflicts we see from the MEDUSA coupling interface developments are with code contained in the GO6 definition and not with the coupling specific code added by the extra GC3 branches. Specifically, we find that dev_r5518_coupling_GSI7_GSI8_landice_bitcomp is a major source of clashes. There is no way to rework the MEDUSA interface to avoid this because this branch, although only part of the GO6 configuration, actually features changes relating to coupled models (for reasons of avoiding code clashes!) particularly w.r.t. the magic numbers used for coupled field indexing.

Discussing with Tim Graham, he's happy for the MEDUSA interface to be added to the GO6 definition (i.e. to dev_r5518_coupling_GSI7_GSI8_landice_bitcomp or a variant of that branch) provided we can disable the MEDUSA code to allow that branch to continue to be compatible with all existing non MEDUSA configurations.

That should be do-able. It seems we'll have to control the relevant code with key_medusa, which is not ideal, unless we can introduce dummy modules containing dummy variable declarations which would allow us to control things within sbccpl.F90, etc purely with run time logicals.

Initially I start on the basis of controlling things with CPP keys.

  • Created new NEMO vn3.6 branch: dev_r5518_GSI7_GSI8_landice_bitcomp_medusa
  • Merged in dev_r5518_coupling_GSI7_GSI8_landice_bitcomp at revision 6659
  • Added MEDUSA-specific coupling interface changes to sbccpl.F90 (revisions 6659-6670 inclusive). NOTE: These are not finalised in any sense since we have no UM/coupler source for the incoming dust, the incoming f_pco2a field is not defined as a 2D field in MEDUSA, merely as a single value scalar and the f_co2flux2d probably has nowhere to go in the UM while bugs remain to be fixed in the UM TRIFFID scheme.
  • So the only field we might reasonably couple and expect to deal with realistically at the moment is the MEDUSA→UM DMS field.


  • Compilation W/ MEDUSA: Using a copy of u-ad903 we have successful compilation. BUT to achieve this I have to take a copy of Julien's dev_r5518_NOC_MEDUSA_Stable and make modifications to ensure that various necessary fields are made publicly available in order to allow them to be referenced in the main NEMO code. The method of doing this is to move the declaration of the necessary fields to the top of the appropriate module definition. This is not necessarily ideal since there are numerous similar variables not required for coupling which remain defined within subroutines in the modules. It may be that an alternative approach is to define the coupling fields within sbccpl or cpl_oasis and then have those referenced by the MEDUSA code. This approach may also allow us to dispense with CPP keys in the main NEMO code but does have the overhead of needing MEDUSA to copy the fields to and from their ultimate target and source arrays. However, for now we stick to CPP keys and our model compiles OK using the following 3 NEMO branches/WC's:
  • /data/local/frrh/nemo_br/dev_r5518_GO6_package - My WC of the GO6 package branch featuring my extra changes from dev_r5518_GSI7_GSI8_landice_bitcomp_medusa
  • branches/UKMO/dev_r5518_GC3_couple_pkg - The GC3 coupling package branch
  • /data/local/frrh/nemo_br/dev_r5518_RH_MEDUSA_Stable_mod - My WC of dev_r5518_NOC_MEDUSA_Stable with additions to make coupling fields available publicly.
  • Compilation W/O MEDUSA: Using a copy of mi-aj246, named mi-aj246_CRUN_fix_MEDcode_Nomed which is just a NEMO_CICE job - I include /data/local/frrh/nemo_br/dev_r5518_GO6_package - My WC of the GO6 package branch featuring my extra changes from dev_r5518_GSI7_GSI8_landice_bitcomp_medusa. This compiles OK but fails in the run with an XIOS error because the iodef.xml file is incompatible with
  • Running W/O coupling:
    • Using suite u-ad903, but with no special changes to the namcouple and no additional changes to the NEMO namsbc_cpl namelist, we have a successful run for 1 day. i.e. This is more or less a standard GC3

set-up with MEDUSA included in the ocean component but no extra Atmos-MEDUSA coupling.

  • Running W/ coupling:
    • We can't activate all the coupling fields because most of them are not ready…. f_pco2a is still a scalar value… it needs to be a 2d array in JP's MEDUSA branch. Dust cannot be coupled as we have no source term in the atmos, pending suitable code developments. MEDUSA CO2 flux is an unknown quantity due to unknown robustness of TRIFFID code. That leaves DMS. We'll try that….
  • Take a WC of our suite… u-ad903_dmsCouple.
  • Modify NEMO namsbc_cpl to add entries for our four new fields with only DMS active.
  • Modify the suite namcouple to add DMS (only).
  • Submit.

  • Well this fails - initially because it was trying to couple fields which I didn't want to exchange. It turns out the best way to manage these is via checking sn_rcv_atm_pco2%cldes, sn_rcv_atm_dust%cldes, sn_snd_bio_dms%cldes and sn_snd_bio_co2%cldes values. We use a value of "medusa" in the suite GUI to enable each field. I use "disable" to switch it off as far as NEMO is concerned, but you could just as well use "tomato" or anything which isn't "medusa".
  • Switching just to using DMS, which I believe is probably the most likely to be viable field, the thing still crashes. It turns out that the DMS field dms_surf2d and f_co2flux2d and f_pco2a are only temporary work arrays in trcbio_medusa. They are allocated, used and deallocated entirely within the life of that routine, so are completely unsuitable for use as coupling fields in the puts and gets.
  • Contacted JP to explain the problems. I could do something along the lines of the following:
      create special coupling arrays which I declare and use in 
      the NEMO coupling interface, say: DMS_out_cpl, CO2Flux_out_cpl, PCO2_in_cpl
      and  Dust_in_cpl. Then the MEDUSA code would have to add things like
       USE sbccpl, ONLY: DMS_out_cpl, CO2Flux_out_cpl
       USE sbccpl, ONLY: PCO2_in_cpl, Dust_in_cpl
       in the appropriate MEDUSA modules and just copy the appropriate values
       to and from these arrays at the correct times. 
    The problem with that is that we don't really know when the coupling exchanges have happened or are about to happen so you could end up copying to and from these arrays on every timestep which is a small cost in terms of the overall model costs, but it's still not a great precedent to set. Well we're net even setting a precedent because sbccpl is full of stuff copying and processing fields on every TS which may or may not ultimately result in an actual PUT operation. So maybe as an approach it's not that bad.

So now I modify oce.F90 to add the four fields indicated above and allocate them via a CPP condition on key_oasis3 and key_medusa. I also introduce new RTL defined in oce.F90 for now : ln_medusa to enable us to test things at run time without littering the code with #ifdefs.

This compiles and fails - we have no values in DMS. The routine where they seem to be available is trcbio_medusa. But we have to intercept values before the arrays are deallocated. It's made extra difficult by the fact that trcbio_medusa basically contains one massive 3 level loop: over k,j,i. It's huge and it's not clear at what stage we should populate our arrays - for now I do it just prior to the deallocation. That's all jolly well but it doesn't solve the issue of where the values come from at startup because this code will not be called until after the first coupling exchange. For now I just initialise DMS_out_cpl and CO2Flux_out_cpl to zero in oce.F90, but that is purely for testing purposes and must not be seen as a solution because it clearly does not cater for restartability.

OK so this fails to run - We need to enable l_oasis_obgc in the suite to ensure arrays are allocated to send/receive coupling variables. Doing that we then gat a reconfiguration failure as we need to include fields 0,197 in the UM start dump. Add this in the "Configure ancils and initialise dump files" section and initialise to zero. That gets through the reconfiguration again. But do we have a source for the incoming DMS beyond the recieved field. I'm not sure we do…. That's probably fine for now, for testing purposes, but where does it need to go ultimately?

Well that now runs OK. (Remember the atmosphere is merely receiving the DMS field, it's not using it anywhere… that's for the code owner to figure out.) Checking the DMS coupling field using EXPOUT we see that at TS=1 we have uniform zero and non zero values evolve at subsequent coupling exchanges. However, looking at the incoming field to the atmos we see a column of an apparent steep gradient in the field below India. This is the point at which the orca grid begins and ends - i.e. this looks like a wrap around issue with the MEDUSA DMS field.

Sure enough, the outgoing DMS field from NEMO has zeros in the first and last columns! This is consistent with the loops in trcbio_medusa.F90 which only operate over columns 2 to 361 and rows 2 to 331. Is that deliberate? One presumes it must be otherwise MEDUSA runs would have bad things going on around the grid edges. If it is all expected, then we probably need to add a halo update before coupling DMS (and CO2 flux).

What I've done for now is to modify a working copy of the MEDUSA branch to add an explicit lbc_lnk to DMS_out_cpl as soon as it's populated in trcbio_medusa. So dev_r5518_GSI7_GSI8_landice_bitcomp_medusa is safe to go into GO6/GC3 now and the only key_medusa refernces are in oce.F90 where they set up ln_medusa which is used at run time and the allocation of coupling fields if key_oasis3 is defined. As a consequence of this approach, no extra keys are needed in sbccpl.F90 as we can control everything by the RTL ln_medusa. However this branch is still work in progress because we need to remove the initialisation of couplig fields in oce.F90 when the MEDUSA branch supplies us with valid initial coupling fields.