Climatology and nudging in ecosystem models

[jklinck]

I have an ecosystem model (a modification from Fennel.h) in which I would like to nudge nutrients to climatology. I am trying to define the climatology in ana_tclima.h with statements like those below. I also include "use mod_biology" in the module section which should give access to the biological variables.

Compile of analytical.f90 fails due to undefined tracer indexes (iNO3_, iNH4_, etc). There is an old example in ana_tclima.h from bestnpz which refers to iFe. Maybe this sort of reference went away with TCLIMATOLOGY. Maybe I need to use the real (integer) tracer indexes (seems a bit clumsy).

So: How do I define climatology for ecosystem variables in ana_tclima.h?

Thanks for the help.

John

CLIMA(ng)%tclm(i,j,k,itemp)=T0(ng)
CLIMA(ng)%tclm(i,j,k,isalt)=S0(ng)
CLIMA(ng)%tclm(i,j,k,iNO3_)=5.0_r8
CLIMA(ng)%tclm(i,j,k,iNH4_)=1.0_r8
CLIMA(ng)%tclm(i,j,k,iPO4_)=1.0_r8

[wilkin]

I suspect you need to update your varinfo.dat and/or bio_Fennel.in. When we introduced the more generic handling of optional nudging for additional tracers those files had to change.

[jklinck]

My "clever" solution to use explicit indexes for the climatology does not work. The climatology arrays don't have the indexes that I expect so the nudging is wrong.

I have gone back to the "proper" way to have a climatology. I could not find any example of a biological model that uses a climatology, although some biological climatology variables are included in varinfo.dat. The bestnpz.h file has some details but uses an old cpp of TCLIMATOLOGY, so I assume that it does not work correctly.

A line in ana_tclima.h defines the oxygen climatology:

CLIMA(ng)%tclm(i,j,k,iOxyg)=50.0_r8

The compile error is
analytical.f90(3452): error #6404: This name does not have a type, and must have an explicit type. [IOXYG]
CLIMA(ng)%tclm(i,j,k,iOxyg)=50.0_r8
-----------------------------------^

I have changed varinfo.dat as follows:

'oxygen' ! Input/Output
'dissolved oxygen concentration'
'millimole_oxygen meter-3' ! [millimole/m3]
'Oxygen, scalar, series'
'ocean_time'
'idTclm(iOxyg)'
'r3dvar'
1.0d0

The bio_var.h include file has lines:
CASE ('idTclm(iOxyg)')
idTclm(iOxyg)=varid
This file is included in mod_ncparam.F. There does not seem to be any extra info needed in mod_biology.F, other than the bio_mod.h file.

Clearly, I need to change something else, but I am not sure what or where. I am working in the COAWST context, but have checked the ROMS
code as well.

Thanks in advance for your help.

John

[jklinck]

The bio input file has the following lines:

LtracerCLM = 12*T
LnudgeTCLM = 12*T

so climatology should be managed and used.

Nudging to climatology does work with the imposed indexes, but the wrong variables are nudged to the wrong values.

John

[kate]

There was a cleaning out of these chunks from varinfo.dat: 'idTclm(iOxyg)' and there is now code in inp_par.F starting with

Code: Select all

!  Set climatology tracers (active and passive) metadata.

for setting the climatology ids, picking new ids at the end of the list.

[jklinck]

I have figured out how nudging to climatology works, by doing some experiments. Here is my understanding, but I am happy to hear another
explanation.

The climatology arrays are created (and indexed) in the order in which climatology is processed.

In my case, I have 16 tracers (T,S, 4 sediment and 10 biology). I want to nudge T and S, so the ocean.in file has these lines.

LtracerCLM == T T ! temperature, salinity, inert
LnudgeTCLM == T T ! temperature, salinity, inert

and ana_tclim.h has the following lines:
CLIMA(ng)%tclm(i,j,k,itemp)=T0(ng)
CLIMA(ng)%tclm(i,j,k,isalt)=S0(ng)

itemp and isalt are set to 1 and 2, respectively.

I also want to nudge one of the bio variables to a value, specifically the dissolved oxygen, which is bio variable 8, with index iOxyg=14. It is the 14th tracer variable (2 + 4 + .

The bio.in file has the following choices for climatology and nudging where the 8th entry (for oxygen) is T:

LtracerCLM == F F F F F F F T F F
LnudgeTCLM == F F F F F F F T F F

and ana_tclim.h has the following line:

CLIMA(ng)%tclm(i,j,k,3)=30.0_r8

So, the third climatology array has the oxygen climatology. The order of stacking for the climatology is controlled by the order of T values in LtracerCLM, which is determined by the order in which bio variables are created. If I want to nudge another bio variable (say, the fifth
one), then the bio.in file has

LtracerCLM == F F F F T F F T F F
LnudgeTCLM == F F F F T F F T F F

and that climatology will be in array 3 and oxygen would move to array 4.

In step3d_t.F, there is a loop over all values in LnudgeTCLM which is a logical array with all tracers, stacked by physics first, then
sediments, then Bio (I think, but I have not tried nudging sediments). Then true entries identify the tracer index of the variables to be nudged, and a counter points to the associated climatology array. This "indirect" indexing scheme avoids having a
large climatology array (all tracer variables) when only a few variables are nudged.

I hope this is helpful.

John

[wilkin]

John,

What you say here all sounds correct, but let me reiterate for others:

If you are nudging to climatological data that are provided in a netcdf file then it is very straightforward. The climatology file need only include the variables you wish to nudge, and they have the same names as the ROMS state variables (as in the initial conditions, or history files).

We made this change to accommodate biological models where often only a few variables are important for nudging (like open ocean nutrients, and oxygen) to avoid having to have ALL bio variables in the climatology file even if they were not going to be used.

As you have noted, the user controls the activation of nudging with the logical flags. The climatology arrays are only allocated for the variables to be nudged (a great economy on memory).

But this means the itrc counter for the climatology arrays differs from the full tracer state.

If you are nudging to climatological data that are provided in a netcdf file this book keeping is done for you.

If you want to set the climatology array values directly in ana_tclima.h then, yes, you need to be careful about the itrc index.

A robust way to do this is to echo the logic in step3d_t.F. Something lie this ...

Code: Select all

      ic=0
!
      DO itrc=1,NT(ng)
!
!  Set compact reduced memory tracer index for nudging coefficients and
!  climatology arrays.
!
        IF (LtracerCLM(itrc,ng).and.LnudgeTCLM(itrc,ng)) THEN
          ic=ic+1
          DO k=1,N(ng)
            DO j=JstrT,JendT
              DO i=IstrT,IendT
                CLIMA(ng)%tclm(i,j,k,ic)=???
              END DO
            END DO
          END DO
        END IF
      END DO
!

This is robust because the order of indices in the logical array LtracerCLM will match the order in the global tracer array.

I'll let you decide if you start the counter at ic=2 so that you don't set temp and salt here.

[pmaccc]

I am having a similar problem to John's. I am trying to read climatology fields (in this case from a NetCDF file instead of specifying them analytically) and use them for nudging to climatology in a strip near the open boundaries. This works fine for the velocities, and for salt and temp, but the run crashes without comment during initialization right at the moment I believe it should be trying to read the first bio variable climatology field. It is reading it from the same .nc file as where it gets salt and temp. It does read the boundary values of the bio variables successfully.

[kate]

I would try recompiling with USE_DEBUG to see if you get more useful messages. You may be running into some memory or file format limit.

[wilkin]

This works fine for the velocities, and for salt and temp, but the run crashes without comment during initialization right at the moment I believe it should be trying to read the first bio variable climatology field. It is reading it from the same .nc file as where it gets salt and temp.

Make sure you have updated your varinfo.dat and ocean.in files, and check that the biology .in file (for which every ecosystem model you are using) has entries for the logical switches LtracerCLM and LnudgeTCLM set to TRUE for the variables you wish to nudge.

varinfo.dat changed to introduce generic handling of climatology arrays (see this ticket)
https://www.myroms.org/projects/src/ticket/627

[pmaccc]

Thanks very much, John and Kate.

Running with USE_DEBUG=on I get errors like:

forrtl: severe (408): fort: (3): Subscript #1 of the array TVAL has value 0 which is less than the lower bound of 1

I have turned on the flags for processing of bio tracer climatology in my bio.in file, as well as the nudging to climatology, and this appears to be handled correctly as reported in the log file.

I could not find any mention of varinfo.dat in ticket 627 (but I did find that ticket very helpful in prodding me to finally make my own python code to create a nudgcoef.nc file, instead of using ana_nudgcoef.h).

Is there some place in varinfo.dat that allows salt and temp climatology to be read? In the climatology section it specifically mentions u and v, but not salt and temp, yet all 4 are processed correctly as climatology. I keep thinking that if I could copy what is done for salt and temp for my bio variables I would be fine. They are all the same size and shape in the _clm file.

[kate]

forrtl: severe (408): fort: (3): Subscript #1 of the array TVAL has value 0 which is less than the lower bound of 1

Does it give a line number and file name?
I had said:

There was a cleaning out of these chunks from varinfo.dat: 'idTclm(iOxyg)' and there is now code in inp_par.F starting with

Code: Select all

!  Set climatology tracers (active and passive) metadata.

for setting the climatology ids, picking new ids at the end of the list.

You don't put tracer climatology into the varinfo.dat file, but you do put in all the variables in the Tvar section:

Code: Select all

'no3'
  'Nitrate'
  'mol/kg'
  'no3, scalar, series'
  'ocean_time'
  'idTvar(ino3)'
  'r3dvar'
  1.0d0

'river_no3'
  'river runoff Nitrate'
  'mol/kg'
  'no3, scalar, series'
  'ocean_time'
  'idRtrc(ino3)'
  'r3dvar'
  1.0d0

For boundary conditions, you can also have these:

Code: Select all

'no3_east'
  'Nitrate eastern boundary condition'
  'mol/kg'
  'no3_east, scalar, series'
  'ocean_time'
  'idTbry(ieast,ino3)'
  'nulvar'
  1.0d0

(one for each side).

Also, it could be failing on the time variable. Best to provide a time attribute in your climatology file for each field.

[pmaccc]

Kate,

Here is the complete debugging error message (actually one of many identical ones, presumably because I am running on 72 cores):

forrtl: severe (408): fort: (3): Subscript #1 of the array TVAL has value 0 which is less than the lower bound of 1

Image PC Routine Line Source
oceanG 0000000001CD51A1 get_cycle_ 224 get_cycle.f90
oceanG 0000000001A1992D inquire_ 277 inquire.f90
oceanG 0000000001985BCE get_3dfld_ 104 get_3dfld.f90
oceanG 0000000001305D64 get_data_ 430 get_data.f90
oceanG 0000000000731B64 initial_ 273 initial.f90
oceanG 0000000000426660 ocean_control_mod 130 ocean_control.f90
oceanG 0000000000425300 MAIN__ 95 master.f90
oceanG 0000000000425126 Unknown Unknown Unknown
libc.so.6 00002AAAABF5ECDD Unknown Unknown Unknown
oceanG 0000000000424FA9 Unknown Unknown Unknown

[kate]

You can look at line 430 of your get_data.f90 to find out what fields are being read. You can look at line 224 of get_cycle.f90 which might actually match what I have:

Code: Select all

223           i=MIN(ntime,Tindex+1)
224           Tend=Tval(i)

ntime is an input to the routine giving the number of records in the file. Could that be zero? Tindex is the record number to start reading based on the model time. I would watch the get_cycle logic in the debugger to see what's what with those variables. Or you can add print statements. Yes, you can end up with 72 printed lines, one from each core. You can debug on fewer cores.

[pmaccc]

Kate, related to your comment about the time variable, in our _clm input file here are what salt (whose climatology is read) and NO3 (whose climatology is not read) look like:

In [8]: print(ds['salt'])
<class 'netCDF4._netCDF4.Variable'>
float64 salt(salt_time, s_rho, eta_rho, xi_rho)
long_name: salinity climatology
units: PSU
unlimited dimensions:
current shape = (2, 40, 381, 174)
filling off

In [9]: print(ds['NO3'])
<class 'netCDF4._netCDF4.Variable'>
float64 NO3(salt_time, s_rho, eta_rho, xi_rho)
long_name: NO3 climatology
units: MicroMolar
unlimited dimensions:
current shape = (2, 40, 381, 174)
filling off

There is also a separate, unused, variable called NO3_time, that has the same values as salt_time.

In varinfo.dat the definition of these variables looks like:

'salt' ! Input/Output
'salinity'
'nondimensional' ! [PSU]
'salinity, scalar, series'
'ocean_time'
'idTvar(isalt)'
'r3dvar'
1.0d0

and

'NO3' ! Input/Output
'nitrate concentration'
'millimole_nitrogen meter-3' ! [millimole/m3]
'NO3, scalar, series'
'ocean_time'
'idTvar(iNO3_)'
'r3dvar'
1.0d0

I wonder if I need to redefine NO3 in the _clm file with NO3_time instead of salt_time...?

[kate]

No, I would add a time attribute to each variable so they look like:

Code: Select all

	double bac(ocean_time, s_rho, eta_rho, xi_rho) ;
		bac:_FillValue = 1.00000001504747e+30 ;
		bac:long_name = "averaged bacteria" ;
		bac:units = "mMol N m-3" ;
		bac:time = "ocean_time" ;
		bac:coordinates = "ocean_time s_rho lon_rho lat_rho" ;
		bac:field = "bacteria, scalar, series" ;

except yours should have "salt_time".

[wilkin]

Is there some place in varinfo.dat that allows salt and temp climatology to be read?

What changed with that ticket was that ROMS now propagates the definitions in varinfo.dat for the tracer state variables (temp, salt, NO3 ...) to corresponding climatology variables without having to define entries for the climatology variables explicitly. So stuff was taken out of varinfo.dat, not added. It's much cleaner this way.

Since you get an error related to TVAL that definitely suggests your time coordinate variable is the problem.

Kate's suggestion is the best approach. If your netcdf file includes an attribute "time" for a variable this tells ROMS to use that, and it overrides any default inherited from varinfo.dat. It is definitely the best way to go (and has the advantage of making your climatology file more self-describing).

The "time" attribute is recommended in the netcdf template Data/ROMS/CDL/clm_ts.cdl.

[pmaccc]

Kate and John,

I added the time attribute to the bio climatology variables in their NetCDF file and this appears to have fixed the problem. Hooray and thanks so much for your help!!

PM