MOM4 Frequently Asked Questions

Here are a number of things you can try:

-- wget, it provides command-line access to the http protocol.

-- curl and snarf are two other command-line, script-friendly applications widely available for *NIX machines, allowing the use of a variety of proxies on various ports.

The mom4-test cases are not well tuned models. Hence, they generally will not result in physically meaningful simulations. Additionally, the models may crash if run for longer than suggested in the distributed run scripts. That is, the test cases are supported ONLY for testing the integrity of the code on various platforms (i.e., does the code compile and run for a few time steps?). PLEASE PLEASE do not take the test cases and write papers based on them. We are not providing models--instead we are providing code.

The reason is that ice model and ocean model use different axis. The ocean model simply just uses the grid axis gridlon_t/gridlat_t of the grid spec file as the diag axis of the output file. The ice model uses the average of the geographical longitude/latitude average along j-index/i-index. When the grid is tripolar grid, these two will have different axis.

Try the following setting, suppose ni, nj is the grid size:

It is probable that you are starting your non-boussinesq run from a restart generated from the Boussinesq. That is the problem. The Boussinesq restart does not have the extra rho_nbt variable saved, so the model bombs when switching to the non-Boussinesq in the middle of your run.

Two options:

There is no fundamental reason that test4 cannot run non-Boussinesq.

At GFDL we can run test4 using 1 processor after increasing stacksize by the command: unlimit stacksize

You may need to change coupler_nml dt_ocean if using the coupler_main to drive the ocean, or dt_ocean in ocean_solo_nml if using the ocean_solo.F90 to drive the ocean. The model time steps generally need to change when changing grid resolution. Also, if you alter the layout of the parallel processors used to run the model, then you may need to change ocean_model_nml layout.

So, in general, the test case namelists will need to be changed both for computational and physical reasons. Please do not assume that all is directly transferrable to your new experiment configuration.

The diag_table sets up the diagnostics to be saved during the model integration. If you wish to add or remove some of the diagnostics, then the diag_table will need to be changed. mom4 provides numerous possible diagnostics, so do not assume those fields that are enabled in the diag_table examples are all that are available.

The field table sets up certain numerical and initial conditions for the tracer fields. Many different advection algorithms are set here, as well as certain idealized initial conditions (e.g., constant tracer). The field_table provided in the mom4 tests may happen to be fine for your purposes. But such assumption is not guaranteed.

The data_table is used to bring in forcing fields from datasets to run the ocean model. The file names and the field names of your forcing fields will generally be different from those in the mom4 tests. Remember that if you wish to run with on_grid=.false. for a particular field, then the dataset must be in spherical grid. The exchange grid algorithm assumes that the source field is on a spherical grid as it interpolates to the mom4 spherical or tripolar grids.

If your model grid is the same as that of a mom4 test, then the answer is "yes".

If your grid is different, then you can use the mom4 forcing data if it is on a spherical grid and you set "on_grid=.false." in the data_table. If the mom4 forcing data is on a tripolar grid, then you cannot use this data since the bilinear horizontal interpolation algorithm can only interpolate from a spherical source grid to a spherical or tripolar target grid. In this case, you must regenerate the forcing by following the prodecure discussed in the preprocessing section of this user guide.

If your model grid is the same as that of the mom4 test, then "yes".

If your grid is different, then "no". You must generate the initial conditions on the model grid prior to running. Note that if one wishes to start the model with constant tracers, then it is possible to do so via the appropriate setting in the field_table. More nontrivial initial conditions must be generated from either one of the idealized_ic in the preprocessing portion of the mom4 distribution, or by remapping realistic fields onto the model grid using the remap_3d functionality (again in the preprocessing code).

All forcing data can be mapped onto the model grid once the grid_spec.nc file has been generated.

a) Idealized 3-D data (initial conditions) can be generated through preprocessing/mom4_prep/idealized_ic

b) Realistic 3-D data (initial condition) can be generated through preprocessing/regrid_3d when the source data is on lon/lat grid.

c) Idealized 2-D forcing data can be generated through preprocessing/mom4_prep/idealized_bc. These idealized forcing data will be placed in the files temp_sfc_restore.nc, salt_sfc_restore.nc, tau.nc and water_flux.nc.

d) Realistic 2-D forcing data can be generated through preprocessing/regrid_2d when the input data is on lon/lat grid.

e) When the input dat is on tripolar grid, preprocessing/regrid can be used to regrid 2-D and 3-D forcing data. Keep in mind that preprocessing/regrid is not fully tested hence may fail on some cases. So try to find input forcing data on regular lon/lat grid. If you have to use preprocessing/regrid and encouter some problems, please send an email to the mom4 email list about your issues and the issues will be addressed as soon as possible.

f) If runoff need to be regridded onto your own grid, preprocessing/runoff_regrid should be used. runoff_regrid can only accept input data sets on regular lon/lat grid. This tool moves the river runoff from a spherical grid onto your model grid, whether spherical or tripolar. It places the river runoff onto the sea point nearest to the land near where the original data was located. Care should be taken to be sure you are satisfied with the results of these steps.

You can use the datasets of mom4 test as input for regridding, but it is not recommended when the grid of the mom4 test is tripolar grid. When the input data is on tripolar grid, you have to use preprocessing/regrid. As mentioned in the previous question, you may encounter some problems when you use this tool. If the original lat-lon datasets are available, those datasets are recommended to use.

You can obtain such data sets on the GFDL NOMADS server http://nomads/nomads/forms/mom4.html.

You need to generate exchage grid when your model is driven by the coupler_main in the directory /coupler. You do not need to generate an exchange grid if you are driving the ocean with ocean_solo.F90 in the directory mom4/drivers.

You can check path_names in the input directory to find out if your test is coupled model or solo model. If coupler/coupler_main.f90 is included in the path_names, your test is a coupled model, otherwise it is a solo model.

You can use preprocessing/generate_grid/make_xgrids to generate exchange grid. The usage of make_xgrids is:

make_xgrids -o ocean_grid.nc -a atmos_grid.nc -l land_grid.nc.

Normally atmos_grid.nc and land_grid.nc are chosen to be the same. In this case, the usage will be:

make_xgrids -o ocean_grid.nc -a atmos_grid.nc.

Assume you already generated the grid for ocean model and we call it ocean_grid.nc. If you decide you want to pick atmosphere grid the same as ocean grid, you can generate the exchange grid by: make_xgrids -o ocean_grid.nc -a ocean_grid.nc

If you want to pick different resolution of atmosphere grid, you can use ocean_grid_generator to generate another grid file and passed the grid file to make_xgrids through option -a .

This message is often the result of the model bombing because of a large time step. Here is what typically happens.

If the above is happening, then here is what I suggest:

1. Examine the initialization portion of the stdout file fms.out to see if mom4 has found any a priori problems with your time steps. Are you on the edge of stability, or safely within the linear constraints for the barotropic, baroclinie, and advection time steps? mom4 performs rudimentary time step analyses that are often sufficient to find time step problems. However, it certainly cannot find all problems before running the model.

2. set diag_freq=1 for all places that you can do so, such as the following places:

<namelist name="ocean_model_nml"> energy_diag_freq=1

These diagnostics may help determine where the NaNs start. The diagnostics may also help to determine if there is simply a time step instability (resolved by reducing time steps) or something wrong with your land/sea mask

Model spin-up is, unfortunately, an art rather than a science. Here are some suggestions:

You need to check if the time axis has the cartesian_axis attribute. If it doesn't, you need to add the the following to your data file:

TIME:cartesian_axis = "T" through ncatted.

You can use ncdump/ncgen but this is harder to do with large files since most editors have limited buffer capacity. ncdump/ncgen are included with the NetCDF library.

ncdump foo.nc > foo.cdl

open foo.cdl with your favorite editor and modify header, in this case add "calendar" attribute to time axis

ncgen -o foo_new.nc < foo.cdl

For daily data, forcing fields are linearly interpolated on each call to "time_interp_external" data to the current model time. For example, assuming daily forcing times are centered at 12Z, for model time Jan 1 0Z, the result is (0.5*Dec31 + 0.5*Jan1). The "test_time_interp_ext" program is defaulted to increment the model time in units of months. You can change this through the namelist. The important thing to verify is that the forcing times as printed in the output are the correct values.

Time encoding in the forcing file, e.g. "days since 0001-01-01" or "hours since 0001-01-01" is converted into the model representation of time appropriate either the "no_leap" or "julian" calendars. We currently do not support the Gregorian calendar in the FMS code. This is the source of your problem. We REQUIRE the "calendar_type" attribute associated with the time axis in order to decode the time data correctly. In your case, the NCEP Reanalysis data is encoded using the Gregorian calendar and the reference date is year 1, leading to a substantial drift over nearly 2000 years. Providing the "julain" calendar_type attribute leads to incorrect time decoding. Until the FMS developers at GFDL (or someone in the community) provides Gregorian calendar support in the FMS "time_manager" module, the only suggestion I can give is to either change the time data/metadata using tools such as the NCO operators, or use an application such as Ferret . Sorry for this inconvenience, but I've been requesting support for the Gregorian calendar for some time (check the GForge site under MOM4 feature requests).

Here's an example Ferret script to save a modified version of the NCEP Reanalysis data with time values which can be handled correctly in the model:

use "ncep_annual_forcing_file.nc"
def ax/t="31-dec-($year) 12:00:00":"02-jan-($year) 
12:00:00":24/units="hours" time_out
define grid/like=var/t=time_out grid_out
let var_mod = var[gt=grid_out@ave]
repeat/l=1:365 save/app/file="var_mod.nc" var_mod
 

The problem can be solved by adding namelist fms_nml in your runscripts.

    &fms_nml
       domains_stack_size = 70656/
 

where 70656 is the number appeared in the error message.

If you want to use fms_io read_data, do NOT pass domain in the list of args, read ALL data onto a temporary array like tmp(1, 90) then in each PE just assign eta(1,jsc:jec)=tmp(1,jsc:jec). Besides, you can also use netcdf routines directly or use mpp_io to do this task.

In ocean_solo.F90, the atmospheric pressure is held inside of the derived type Ice_ocean_boundary%p

This field is read from the data_table entry just as you suggest. The atmospheric pressure is then read into the ocean model field "patm" inside ocean_sbc.F90 subroutine using the code

   if(use_waterflux) then
      var = grav*rho0*max_ice_thickness
      do j = jsc_bnd,jec_bnd
         do i = isc_bnd,iec_bnd
            ii = i + i_shift
            jj = j + j_shift
            patm(ii,jj) = min(Ice_ocean_boundary%p(i,j),var)
         enddo
      enddo
      call mpp_update_domains(patm(:,:), Dom%domain2d)
   endif

We have this min function as a cludge to keep the top model grid cell from being lost in cases where the ice is very thick.

Note the need to run with use_waterflux=.true. to get a nontrivial patm into mom4. This constraint is necessitated when using an ice model in order to get water balances correct for the combined system ocean+ice. If you are using use_waterflux=.false., then you can modify the code as need be. But I recommend against this change.

So bottomline, your data_table is correct. But to get a nontrivial patm into mom4p0c, you need to have use_waterflux=.true.

During the course of model integration user may want to replace a model field with values from a data file or even with a const. Data_table lists all fields you want to override and the source of data used for override. Data_table is also a convenient way to initialize model. More details are available in shared/data_override.

Diag_table is for writing model's results. The table lists all field names, file names, output frequency and other information. More details are available in shared/diag_manager.

Let's call the time at the top of diag_table T1 and the model's initial time T2. The time at the top of the diag_table (T1) is used as starting time of time axis in all output files that have time axis. When you open any output file and look for variable Time, you will see:

      Time:units = "hours since T1"
    

(assuming unit is hours).

T2 is the initial time of the model, all fields cannot be written to output files sooner than T2. Since time can not be negative you should have T2 >= T1. In practice, to keep values in time axis from being exceedingly large, you'd better have T2=T1.

It is possible that in the diag_table you have initial time = Dec. 31 19xx and the saving frequency is MONTHS. The day is fixed (31) and the month is incremented by 1. You will end up having Feb. 31 which is illegal, hence the error message.

The xland portion in the field table specifies the rate of flow between two points, one inside and one outside of a land-locked region. That is, we mix between (i,j) points (ixland_1,jxland_1) and (ixland_2,jxland_2).

In the example below, there are two sets of points chosen: set A and set B.

"xland","Gibraltar","ixland_1=274,ixland_2=276,jxland_1=146,jxland_2=146,kxland_1=1,kxland_2=28,vxland=0.55e6"

set A mixes between (i,j) point (274,146) and (276,146) over k-levels 1,28 with rate 0.55e6 m^3/sec

"xland","Gibraltar","ixland_1=274,ixland_2=276,jxland_1=147,jxland_2=147,kxland_1=1,kxland_2=28,vxland=0.55e6"

set B mixes between (i,j) point (274,147) and (276,147) over k-levels 1,28 with rate 0.55e6 m^3/sec

Note that set B is one latitude row north of set A. We recommend mixing in coupled sets, such as those chosen here, in order to minimize the splitting between grid cells that can occur on the B-grid. That is, it is important to couple the black and white squares on the "checkerboard". This point is mentioned in Section 13.8 of the mom4 technical guide.

All tracers are affected similarly. And since k=1 is chosen, volume is mixed as well as tracer.

At GFDL, we try our best to make the scripts independent from platforms as much as possible. In theory, all platform specifics are contained in Make templates (mkmf.template.platform). When users set platform (sgi, ibm, ifc ...) the scripts will pick the right template.

In reality, users encounter various problems related to implementation on different platforms. Users may need to modify the source code and/or scripts. Please refer to MOM4 on Various Computer Platforms section to find if your problems are listed there. If not, please report your problems to us and we try our best to address your problems.

In the runscript you see:

      set days = 1
      set months = 0
    

You can change days and months to any number you want. However, note that the tests are provided for checking ports of the code to different machines, and overall numerical integrity of the code. The tests are not guaranteed to run for an indefinite period of time (i.e., they are not `` tuned models''). A test case therefore may bomb if changing the runlength to something longer.

In the runscript you see:

      set npes = 15
    

In this example, there are 15 processors, you can set npes to any number your system can afford. However, You need to pay attention in case of STATIC option.

Based on our experience with SGI Origins at GFDL, The compiler is performing an excessive number of array copies if all the arrays are dynamically allocated. The solution is to allocate arrays at compile time.

Requirements of static option: each processor should have the same number of grid points as all other processor. For example, if your global grid is 100x100 you can not have 13 processors but you can have 2, 4, 8, 10 ... processors.

For more details see the script for static option (test1/run_mom4_test1_static).

In the runscript look for &ocean_solo_nml (or &ocean_coupler_nml if you run coupled model). You will see:

      dt_ocean = 3600
    

This is time step of ocean model measured in seconds. Note that a longer time step may violate CFL conditions and the model may bomb.

"dt_ocean" corresponds to the ocean tracer time step. The other model time steps are set relative to one another via the correspondence

       dt_ocean = dtts
                = baroclinic_split * dtuv
                = surface_height_split * dteta
                = barotropic_split * baroclinic_split * dtfs
    

with

       dtts  = tracer time step
       dtuv  = baroclinic momentum time step
       dteta = "big leap-frog" time step (the time which updates the surface height)
       dtfs  = barotropic momentum time step
    

Output files are compressed in tar format. First, you need to uncompress the output files. The command for extracting a tar file is:

      tar -xvf filename.tar
    

As a result, you will get a number of files in NetCDF format. Some common tools working with NetCDF are ncview and ferret.

Basically, you need to use restart data instead of initial condition data. If you are using the runscripts provided by the mom4p0 release. Do the following:

      Remove the following section from the runscript:

         cp $expdir/preprocessing/ocean_tracers_ic.tar .
         tar -xvf ocean_tracers_ic.tar
         rm -f ocean_tracers_ic.tar

      and add the following in the same place:

         cp restart_file.nc.tar . (including the path)
         tar -xvf restart_file.nc.tar
         rm -f restart_file.nc.tar
    

It is good idea for each test case to check if:

    - model results reproduce across different sets of processors
    - model results reproduce across restart data
    

In both cases users can rely on checksums that are printed at the end of each run (output file fms.out). For example a run using 1 processor and another run using 10 processors should have the same checksums. NetCDF tool ncdiff can also be used to view differences of two .nc files.

MOM4 has plenty of material written to stdout (which is what is written to fms.out). The user is strongly recommended to become familiar with the contents of fms.out, especially when building a new model experiment. This file contains numerous notes, warnings, error messages, and numerical diagnostics. Our experience is that most of the problems running MOM4 can be resolved by reading through the fms.out file.

When running the model on many processors each processor will have a separate output file (like file.nc.0000, file.nc.0001, ...). mppnccombine will put all these separate pieces together so that we have a file corresponding to the whole grid.

Here is the combining command: $MPPNCCOMBINE file.nc where file.nc is the unified output file, the input is not specified in the command line because it is assumed to be identical to the output plus 0000, 0001, .... By the way running the combining command without any argument will give help on how to use this command

The warning takes place only when truncate_eta=.true. in ocean_freesurf_nml. This is intended for spinups where the surface height may become extremely large and result in negative thickness for the top model grid cell. Truncating eta to a small value will NOT conserve ocean properties. Hence, this option is not recommended for experiments where you are interested in the results.

In ocean_core/ocean_sbc.F90, search for "lprec". You will see that lprec is part of Ice_ocean_boundary derived type. The assumed units of lprec are mass units, so there is a division in the code by rho0 to bring it to m/s. Here are the comments from ocean_sbc.F90

! freshwater flux in kg/(m^2 sec) from rivers and calving land glaciers ! rho0r converts mass flux to (m/s)

It is possible that the output frequency (defined in diag_table) is larger than the model's run length. If output frequency is a month and the model has run only 10 days there will be zeroes in model output. You need to change either output frequency or run length.

The namelist "use_xlandmix" refers to the algorithm "mom4/ocean_param/sources/xlandmix/ocean_xlandmix.f90" which handles isolated basins such as the Mediterranean, Red Sea, etc. whose relatively narrow passages are not resolved by the model resolution. In these cases, an ASCII formatted table containing mixing rates and a depth range between selected points in grid ccordinates can be supplied to achieve the desired exchange of tracers. You can refer to the html file in the above directory for more information. If you are not interested in achieving this effect or your model does not contain isolated basins, then you can safely remove the table from your runscript. Otherwise, you should consider examining the table entries and enabling this option.

When you generate topography, there is an option

fill_first_row

My guess is that you have this set .true. in your preprocessing topography generation. In this case, the zeros that you are reading at j=1 mean this is land. The reason that we have this option is due to technical details related to coupling to our ice model. If you are not coupling to the GFDL sea ice model, then fill_first_row=.false. is fine.

Here are some things to consider:

-- are you indeed using the constant vertical diffusion module, or perhaps are you using ppvmix or kppvmix but do not realize this to be the case?

-- do you see vertically unstratified water in your solution?

-- why would you expect to get a zero overturning without convection? There are still downwelling and upwelling near boundaries that may contribute to vertical motion.

We always have frazil enabled (see mom4/ocean_core/ocean_tracer.F90) in order to allow the diagnostic tracer type "T_diag" to be nontrivial, even with idealized simulations where frazil is not used. It is a bit unclean, but there are trivial computations added which should be of minimal consequence to the cpu cost.

The mom4 default has no latitudinal dependence, since that is the traditional implementation of the scheme according to Bryan and Lewis (1979). We have found some utility in our coupled modelling with increased vertical mixing in the higher latitudes, especially the NAtl. Hence the new flexibility added to mom4.

Please refer to section 5.6 of the online User's manual for more details. In general, the test cases need some kind of "tuning" to make them more scalable. The modification may involve diag_freq, diag_table, clock_grain in fms_nml. Also, one should not count times that are spent in initialization and termination of a run since these may involve heavy I/O activities

To create the overturning streamfunction, you must create one from the following diagnostics added to the diag_table:

tx_trans
ty_trans
tx_trans_rho
ty_trans_rho
ty_trans_gm
tx_trans_gm
 

These fields have the grid factors built in, and they have dimensions of Sv. The code computing these fields is in ocean_adv_vel_diag.F90 and ocean_neutral_physics.F90.

To get global overturning in (j,k) space using Ferret, say

let moc = ty_trans[i=@sum,k=@rsum] + ty_trans_gm[i=@sum]
 

Note the difference in sums. There is no indefinite sum needed for ty_trans_gm, since the vertical integral is already performed analytically within mom4. If you wish to see what I mean, then look at equation (9.97) in Stephen Griffies' book.

If you have a basin mask defined (mask=1 in region of interest, 0 elsewhere), then you can construct by-basin streamfunctions within Ferret using the following commands.

let trans_atl_euler = ty_trans*tmask_atl
let trans_atl_gm    = ty_trans_gm*tmask_atl
let moc_atl = trans_atl_euler[i=@sum,k=@rsum] + trans_atl_gm[i=@sum]
 

You can construct overturning mapped to potential density surfaces by using ty_trans_rho. Note that this diagnostic only incorporates the resolved velocity field. We do not have the GM velocity mapped to density surfaces.

You can obtain such information by comparing global and compute domain. You can get compute and global domain by pseudocode like this:

  call mpp_get_compute_domain(domain, isc, iec, jsc, jec)
  call mpp_get_global_domain(domain, isg, ieg, jsg, jeg)

  if(isc == isg) westmost
  if(iec == ieg) eastmost
  if(jsc == jsg) southmost
  if(jec == jeg) northmost
 

use: LDFLAGS= -brename:.flush,.flush_ -brename:.memuse,.memuse_ in your makefile or in the mkmf template.

There's a typo in your compiler flag, it should contain the following: LDFLAGS= -brename:.memuse,.memuse_

This is probably because memuse.c is compiled with 32 bit option while other files are compiled with 64 bit. You may need to add CFLAGS=-64 to your SGI template. This will ensure that memuse.o will be a 64-bit object.

The problem is correctly diagnosed. A quick solution is to simply remove the "bounds" attribute from the dataset using the NCO utility "ncatted". In the future, MOM4 and Ferret developers will cooperate to solve this problem.