From: Helge Goessling
Sent: 10 March 2021 12:07
Subject: Re: Spectral nuding / relaxation code in OpenIFS
Dear Markus and Glenn,
I am coming back to the question of spectral nudging with OpenIFS. In the meantime we have had some discussions with Etienne Tourigny and Vladimir Lapin at BSC (they do grid-point nudging, so far with OpenIFS-cy40 and an earlier IFS cycle in ECEarth3, if I’m not mistaken), and with Lauriane Batté at Météo France where they do spectral nudging, but in ARPEGE. We still have some questions remaining, and maybe you can help us:
1) Do you know anyone that has been doing spectral nudging with OpenIFS (ideally with cy43)?
2) As mentioned in our earlier exchanges (see below), we have identified the files ifs/adiab/spchor.F90, ifs/module/yomnud.F90, ifs/setup/sunud.F90, and the corresponding namelist namnud to be central to the spectral nudging (they are the counterparts to corresponding files that exist for the grid-point nudging). The same set of files is also used in ARPEGE for the spectral nudging. Thus we assume that in principle the spectral nudging should be functioning in OpenIFS. Would you agree, or are there parts of that code that might not have been updated to account for other changes in the code, or similar, so that some extra coding would be required to make this working (again)?
3) Do the spectral nudging and the grid-point nudging require the exact same nudging files (that is, primarily spectral fields in spectral space, and primarily grid-point fields in grid-point space)?
4) Regarding the preparation of the nudging files (from ERA5), the ECMWF confluence page explains a technique involving PrepIFS. In contrast, both the BSC and MF colleagues do this in a different way that involves running the model in a special mode for individual timesteps. To make things even more diverse, we are considering trying to achieve this with CDO, which includes the operator remapeta to interpolate vertically from one set of hybrid model levels to another set of hybrid model levels. That is working fine with ECHAM, so we hope this can be transferred to OpenIFS. We would be very curious to hear from you what you think about these different options.
Any help would be great!!
Cheers,
Helge
Am 27.07.2020 um 06:33 schrieb Thomas Jung:
Hi all,
yes, that looks very much like the code I was having in mind. I used this as a blueprint for adding grid point nudging or relaxation to the IFS many years ago... ;-)
Best wishes,
Thomas
On 26.07.20 17:25, Glenn Carver wrote:
Hi Jan,
That looks like it. Whether it still works is a question though. Marcus might know better than me.
I think the comments are correct, there is likely a smooth transition between nudged and unnudged wavenumbers.
Cheers, Glenn
From: Jan Streffing
Sent: 26 July 2020 12:41
Subject: Re: Spectral nuding / relaxation code in OpenIFS
Hello toghether,
I went looking for a wavelength limiter for spectral mapping and I found in the code:
ifs/module/yomnud.F90:! NSPNU1 : wave number below which full nudging is applied
ifs/module/yomnud.F90:! NSPNU2 : wave number beyond which no nudging is applied
That seems almost right, although from the wording of the comment they would apprear to do the same thing. I suspect one of the descriptions might be misleading and these are actually what Thomas is looking for. These variables were introduced by MeteoFrance in cycle 32.
These two finally make their way into
ifs/adiab/spchor.F90
Where they are applied as following: (example from Vorticity nudging)
!* 3 NUDGING.
! --------
IF (LNUDG) THEN
! 3.1 GMV VARIABLES:
! * Vorticity:
IF(LNUDVO)THEN
IF(XNUDVO <= 0) THEN
DO JSP=KSTA,KEND
ZUSNUD2(JSP)=XNUDVO
ZUSNUD1(JSP)=1.0_JPRB
ENDDO
ELSE
DO JSP=KSTA,KEND
ZUSNUD2(JSP)=XNUDVO*XWNUDG*&
&MIN(1.0_JPRB,MAX(0.0_JPRB,1.0_JPRB*(NSPNU2-YDLAP%NVALUE(JSP))/REAL(NSPNU2-NSPNU1)))
ZUSNUD1(JSP)=1.0_JPRB/(1.0_JPRB+ZUSNUD2(JSP))
ENDDO
ENDIF
DO JLEV=1,NFLEVL
ILEV=MYLEVS(JLEV)
IF (ILEV >= NTNUDG) THEN
ZREF(:)=0.0_JPRB
DO JSTEP=1,NFNUDG
DO JSP=KSTA,KEND
ZREF(JSP)=ZREF(JSP)+TNUDVO(JSP,JLEV,JSTEP)*XPNUDG(JSTEP)
ENDDO
ENDDO
DO JSP=KSTA,KEND
PSPVOR(JLEV,JSP)=(PSPVOR(JLEV,JSP)+ZUSNUD2(JSP)*ZREF(JSP))*ZUSNUD1(JSP)
ENDDO
ENDIF
ENDDO
ENDIF
Cheers, Jan
On 26.07.20 12:39, Thomas Jung wrote:
Dear All,
it is great to see that there is so much interest in the nudging code. From what I can see, this formulation does nudging / relaxation in grid point space. Having this option is great, since it allows to localize nudging regionally. However, there are also good reasons for doing nudging in spectral space, since this will allow to constrain only the longer waves (e.g. the jet stream). From what I can see this is currently not implemented (in this case I would expect to see something in the namelist regarding total wavenumbers). To my understanding, however, there is a spectral nudging version in the IFS (or older cycles) that has been put in by MeteoFrance. In fact, given that OpenIFS is s spectral model, being able to exploit this capability diagnostically would be a major asset.
Best wishes,
Thomas
On 26.07.20 12:30, Glenn Carver wrote:
Hi Jan,
Just to add that no nudging code was removed between IFS and OpenIFS. The changes to the namelist that Marcus mentioned are in the next release of 43r3 which will be in a couple of weeks.
Glenn
From: Jan Streffing
Sent: 25 July 2020 12:45
Subject: Re: Spectral nuding / relaxation code in OpenIFS
Hello Marcus,
thank you for sharing you knowledge. I have included the interested parties from AWI in cc. Maybe Helge can elaborate on his ideas regarding different nudging schemes and whether what OpenIFS 43 inherited from IFS 43 is what he would like to use.
Best regards, Jan
16 Comments
Marcus Koehler
Dear Helge and colleagues,
As far as I am aware there are indeed two separate nudging schemes in the OpenIFS model: NAMRLX to control nudging in grid point space (reading both spectral and grid point relaxation files) and NAMNUD which controls nudging by constraining wave numbers (as I believe). I don't know about the origin of these two methods and I have so far only used the NAMRLX nudging code.
The Confluence page on How to use OpenIFS in nudged configuration refers to nudging in grid point space with NAMRLX.
With regard to your original questions at the beginning of the above email trail:
1) Unfortunately I am not aware of anyone who uses OpenIFS yet with spectral nudging. I will invite colleagues from the seasonal forecasting to look at this thread to see if they have used it with the IFS.
2)+3) While I haven't used the NAMNUD code myself as yet, I don't know of any reason why this code should not be working in OpenIFS 43r3. As Unknown User (nagc) wrote earlier (in the email thread above) none of the original code was removed during the conversion from IFS to OpenIFS. I have so far identified about 10 source code files that are linked to this nudging scheme:
src/ifs/climate/updtim.F90
Unfortunately I haven't had the chance yet to familiarise myself with how this works and what format the relaxation files need to have. If you know anyone at MF or from the ARPEGE community who has worked with this then I am sure they would be able to advise about the relaxation file format. I noticed that the NAMNUD namelist provides many more switches for nudging of surface boundary parameters (which have no spectral representation) compared to NAMRLX.
4) At ECMWF we can produce nudging files for the NAMRLX nudging code, however for the time being this is only based on ERA-Interim at T255 and on ERA5 at T639. The suite that we use for this does not (yet) include the interpolation to a specified model grid resolution.
Meanwhile we are working on a web-based facility to produce initial data and nudging files, however this will only become available later this year (after summer 2021).
I am aware that the colleagues at BSC have developed their own set of scripts to produce nudging files for OpenIFS 40r1 and as far as I know they use OpenIFS with fullpos to interpolate to the desired output grid. I believe they are working now on a version for OpenIFS 43r3 because cycle 40r1 nudging files may not necessarily work with cycle 43r3 due to the changes in the land sea mask and other changes.
I am not sure about using CDO. I believe this may probably work fine with vertical levels however I can't really say anything about the horizontal grid (especially as OpenIFS 43r3 can also use the cubic octahedral grid).
Sorry if this is only of limited help.
Best regards,
Marcus
Tim Stockdale
Starting in Cycle 43R1 of the IFS, NAMRLX includes NRLXSMAX, which limits the nudging of spectral fields up to a certain spectral wavenumber. Unlike YOMNUD, there is a sharp cut-off. I also have a more recent branch which has an option for zonal mean only nudging. It is relatively easy to tweak the source code to control which wavenumbers of spectal fields are nudged using the NAMRLX/grid point implementation of nudging, so there is no need to use the NAMNUD code if the only requirement is to apply spectral truncation. (I have no experience of NAMNUD< so cannot comment on what else it might do, e.g. with surface fields. I would want to be sure, though, that it was acting consistently with the IFS physics and timestepping, and not just ARPEGE).
In RD we have poorly documented but effective interpolation code, that is both accurate and vastly more efficient than using IFS configurations (presumably in FULPOS mode) to do the interpolation. From its documentation, the CDO remapeta looks as if it is designed to do a sensible and well designed vertical interpolation. I don't know how efficient it would be for processing large data volumes, and it looks as if it would need to be configured as part of a workflow with first a horizontal interpolation, and making sure all appropriate fields are available and in the right format. (TBH, the RD code also has this requirement, which is one reason the whole process is a bit complicated). The only thing that the RD code has that remapeta seems not to is the option for simultaneously perfectly conserving and highly smooth interpolation, which I needed for some specific ozone work but is not at all needed for general use. If an effective workflow can be set up using CDO, the right orography files are specified, and the run-time is acceptable, I would think remapeta looks a good way to go. The results would have to be checked, of course, that everything works OK, particularly where there are large changes in orography.
Helge Goessling
Thanks, Marcus Koehler and Tim Stockdale , this is extremely helpful!
In ECHAM we are also using a sharp cut-off in terms of wavenumber, and as we basically just want to transfer the approach used in ECHAM to OpenIFS, it seems like there's no need to further try using NAMNUD, but do everything simply with the well-tested NAMRLX and specifying NRLXSMAX. I must however admit that I am a bit confused how the spectral truncation goes together with constraining the nudging to certain geographical regions, which to my understanding can be done with NAMRLX not only for grid-point variable but also for the (originally) spectral variables. Tim Stockdale, could you comment on this?
Regarding remapeta in CDO, indeed horizontal transformation (sp2gp and gp2sp) and interpolation need to be done before and afterwards, which are also all coded in other CDO operators, and as far as I know, in the newest CDO versions they also work with the TCo grids. So, also considering your thoughts, I think I remain optimism that this might turn out to be a good way to prepare OpenIFS nudging files.
Tim Stockdale
Yes, it is strange to think of selecting certain wavenumbers and imposing a geographical region at the same time! The way it works is that first the nudging is applied in spectral space, to some or all of the spectral coefficients. Then the nudged fields are transformed to grid point space, where the original un-nudged fields are also available. The final field is then a combination of the nudged field within a geographical region, and the original unnudged field elsewhere. There is some smoothing applied across a transition zone, and namelist variables control the width of the transition. For high spectral wavenumbers, the result is as you would expect - nudging within the region, no nudging outside. But of course at low wavenumbers, the use of a geographical region will alter the spectrum of the applied nudging. I've never used this myself, but I think it kind of makes sense, if you view the geographical filter as having priority. For example, one could nudge the synoptic but not the mesoscale over the Pacific basin, for example. Anyway, all the best with CDO!
Vladimir Lapin
Many thanks to Marcus Koehler and Tim Stockdale for taking the time to help this newly-formed community of people interested in various nudging capabilities of OIFS.
Helge Goessling, if you do not mind the intrusion, here is my understanding of how spectral truncation and regional masks can be combined in the code. In relaxgp.F90 there is a step where a mask is constructed (ZMASK variable) and then it is applied to all variables, including spectral ones. For SH variables there are grid-space buffer arrays PRLXVO, PRLXDI, PRLXTE, PRLXLP. And grid-point relaxation is applied to those variables after a CALL INV_TRANS(...) which is the interface routine for the inverse spectral transform.
So, doing a truncation before incurring CALL INV_TRANS would effectively add spectral nudging on top grid-point relaxation with just a few lines of code. However, combining willy-nilly sharp regional masks and sharp spectral truncation can surely give rise to undesired artefacts (similar to Gibbs phenomenon).
Etienne Tourigny
One comment about interpolating the gridpoint fields with cdo, I am not sure it is possible to write output in the reduced gaussian grid (which is required by ifs)
So there has to be first a conversion to regular gaussian, but then how would you convert it to reduced gaussian so it can be ingested by ifs?
Etienne Tourigny
Thanks to all for this very interesting discussion!
Replying to myself... it seems you are interested only in the fields stored in spectral files (i.e. vo d t) but not in those stored in gridpoint files (i.e. q ciwc clwc cc o3), so this issue might not be relevant for you.
I assume that if you only apply nudging to the spectral fields the gridpoint files are not even required!
In any case, I think that emos_tool could be used to convert the interpolated gridpoint files in regular gaussian grib format to reduced gaussian grib with the --reduced argument. And you need to do some metadata manipulation of any file processed by cdo (such as setting the gridType with grib_set).
At BSC we have implemented a solution inspired by the workflow of prepIFS experiments (which requires using ifs in a special mode, calling fullpos) to generate Initial Conditions from ECMWF reanalyses, for initializing seasonal predictions.
We use the same method to generate the nudging reference files at 6h timesteps. This is somewhat of an overkill, especially since it requires downloading far more data than the nudging reference files (you need everything to run ifs at every 6h), and requires more time and processing power, as Tim Stockdale says.
So we would be quite interested in testing out a workflow which would only require cdo, eccodes tools (grib_set, etc) and probably emos_tool, especially if it can be made to work for the gridpoint files as well!
Helge Goessling
Etienne Tourigny , I just got the information from Uwe Schulzweida that since CDO-1.9.9 (released November 2020) it is indeed possible to remap to the reduced Gaussian grids. I successfully tried this with both remapnn (nearest neighbour) and remapbil (bilinear), of which I think the latter should be preferred to avoid aliasing. Below is an example that works, using the original nudging grid-point input file available here: http://download.ecmwf.int/test-data/openifs/nudging/examples/ . The first command does linear interpolation to the full grid (replacing 'setgridtype,regular' by 'setgridtype,regularnn' would do nearest-neighhbour interpolation, but again that's probably worse) and converts to NetCDF. The second command interpolates and converts back, again with linear interpolation. The third command summarizes the result. Comparing that to the 'sinfo' output for the original file suggests that the files are structurally very similar; the only difference visible from the summary is that the Dtype of all variables is now 'P16', whereas it was 'P8' for one of the five original variables. There might of course also be other differences in the meta data etc. not revealed by 'sinfo', so some adjustments might need to be done to make this again usable by OpenIFS. We can't test that as long as we have the nudging not yet up and running, so I am curious whether that works at your end
.
> cdo -f nc setgridtype,regular rlxmlgg201308251200 rlxmlgg201308251200_regular.nc cdo setgridtype: Processed 26651400 values from 5 variables over 1 timestep [3.93s 44MB]. > cdo -f grb1 remapbil,rlxmlgg201308251200 rlxmlgg201308251200_regular.nc rlxmlgg201308251200_regular2reduced.grb1 cdo remapbil: Bilinear weights from gaussian (512x256) to gaussian_reduced (88838) grid cdo remapbil: Processed 39321600 values from 5 variables over 1 timestep [0.74s 92MB]. > cdo -sinfo rlxmlgg201308251200_regular2reduced.grb1 File format : GRIB -1 : Institut Source T Steptype Levels Num Points Num Dtype : Parameter ID 1 : ECMWF unknown v instant 60 1 88838 1 P16 : 133.128 2 : ECMWF unknown v instant 60 1 88838 1 P16 : 247.128 3 : ECMWF unknown v instant 60 1 88838 1 P16 : 246.128 4 : ECMWF unknown v instant 60 1 88838 1 P16 : 248.128 5 : ECMWF unknown v instant 60 1 88838 1 P16 : 203.128 Grid coordinates : 1 : gaussian_reduced : points=88838 nlat=256 N128 lat : 89.46282 to -89.46282 degrees_north Vertical coordinates : 1 : hybrid : levels=60 lev : 1 to 60 by 1 level available : vct Time coordinate : 1 step RefTime = 2013-08-25 12:00:00 Units = hours Calendar = proleptic_gregorian YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss YYYY-MM-DD hh:mm:ss 2013-08-25 12:00:00 cdo sinfo: Processed 5 variables over 1 timestep [0.04s 32MB].cdoEtienne Tourigny
Very interesting that cdo can now write to reduced gaussian grid, does that mean that you can use data in reduced gaussian grid in any of the operators? I guess this is not possible if they require the data to be in netcdf format and regular lat/lon grid (which may be the case for the remapeta command).
You can change the grib output precision with `-b P8` but also making sure the netcdf precision is compatible.
Also, a `cdo diff` would be best to make sure there are no large errors in the conversion.
Finally, `grib_ls` and `grib_dump` can show the grib metadata, you will have a difference in the gridType which is easily changed with `grib_set`.
Helge Goessling
Good point, it might well be possible to use remapeta and other CDO operators directly with the reduced-Gaussian data, and probably also the intermediate conversion to netcdf might not be required. I guess either you or we will find that out soon...
And thanks for the other hints!
Helge Goessling
Dear OpenIFS experts,
We are looking for a way to output the nudging increments, so that one can diagnose what the nudging is actually doing. We hope that this way one might be able to better understand what is happening, e.g., in which direction the model tends to "drift", and to optimize the nudging parameters (truncation, timescale, etc.; configurations with the smallest average absolute increments should be preferable?).
We have not found any clues (also in the code) for outputting nudging increments. Does someone know whether that's already possible - and if not, how hard would it be to implement such additional output?
I can imagine that other groups using OpenIFS with nudging might also be interested in such output.
Thanks!
Helge
Tim Stockdale
I don't think there is a way to output the accumulated increments, and it might be a significant job to implement the code to do this properly and cleanly.
However, as long as you are happy with an estimate (which in many cases might be very accurate), then all is not lost. All you need to do is write out the instantaneous values of the nudged field you are nudging at the appropriate frequency (say every 6h). Then take the corresponding values of the field you are nudging to. Take the difference, and use the nudging strength to diagnose what the increment must have been at that timestep. Make sure you work on the original grid. Average all the 6h values to get the estimated mean increment. Work on the original model grid, and use the correct spectral truncations if needed.
Limitations: This will not work well for fields with e.g. a strong diurnal cycle or other form of variability that means 6h sampling gives the wrong mean. It will not work well if you need the detailed time evolution over a short period of time - I use it to get monthly mean values and climatologies of the nudging. Finally, I typically use a 12h nudging, while nudging to data which is linearly interpolated between 6h instantaneous values. Thus my nudging will be relatively smooth over time, and values at the "instantaneous" points will be representative of the time averaged increments. If you use very strong nudging (I have seen 1h used), then it is possible that the nudging in the first timestep after each 6h timepoint may not be so representative of the time average.
For typical upper air fields and 12h nudging, though, this method works well.
You are not the only person to ask about this, and maybe we should see if it is possible to add to the code. The easiest would be if we can add it at a higher calling level and use the existing framework in callpar. I'm not sure this is possible, though (it might get combined with other things that can't be separated). If someone (who knows this better than me) adds code to pass all the necessary memory-distributed arrays down to the relevant routine, we would also need to make sure that this didn't affect the efficiency of the code in the 99.9999% of the cases when nudging is not taking place.
Helge Goessling
Great, thanks a lot for this super fast and helpful response, Tim Stockdale!
I guess we will do this exactly the way you are suggesting it, keeping in mind the sampling limitations. But I agree that this should be fairly accurate in most cases.
If someone would at some point extend the code so that the increments could be output in a clean way, that would of course still be very interesting, to get even more precise information when the diurnal cycle plays a role etc.. Speaking of which, in that case I could imagine that it might be useful to also store accumulated absolute increments (in addition to the normal accumulated increments), so that negative and positive increments e.g. over parts of a diurnal cycle do not cancel each other, if that's possible ... but that would probably just be the icing on the cake.
Vladimir Lapin
Just FYI, I have tried running EC-Earth3 (IFS cycle 36) with ARPEGE's spectral nudging activated (NAMNUD). It did not work and I got
from a dummy routine (ifs-36r4/src/dummy/source/faitou.F90). I think it was called inside src/ifs/climate/updnud.F90. I have no clue why most of the parts of spectral nudging code from ARPEGE are available IFS, but faitou.F90 (and maybe some other routines) were replaced with dummy routines. The situation is the same in OpenIFS-cy43.
But the bottom line is that spectral nudging using NAMNUD is not functional.
Unknown User (gdcarver113@outlook.com)
When EC-Earth3 used IFS alot of the code base that was considered unnecessary was removed. For OpenIFS, code that was not used by ECMWF and therefore unsupported was also removed, this included alot of the MeteoFrance code. Bear in mind that there is also alot of legacy code in IFS that is no longer maintained and unlikely to work, that's also why code was removed from OpenIFS.
From the conversation above, it appears NAMRLX is the way to go for relaxing.
Vladimir Lapin
Thanks for the comment, Ryan. Yes, NAMRLX is the way to go. I decided to try NAMNUD only because the current workhorse version of EC-Earth has IFS cycle 36 which does not allow spectral truncation and applies nudging uniformly across all wavenumbers. And I wanted to add this info to the conversation.