cvmix / cvmix-src Goto Github PK

We need to pass either N_iface or Nsqr_iface to the routine but we want to do computations at cell center. As discussed on the weekly telecon, we want to do all averaging to N_iface the latter is not appropriate for grids with variable vertical resolutions... this means changing

N_cntr(kt)=sqrt((max(Nsqr_iface(kt),cvmix_zero) +                 &
                 max(Nsqr_iface(kt+1),cvmix_zero)) *              &
                 0.5_cvmix_r8)

to

N_cntr(kt)=(sqrt(max(Nsqr_iface(kt),cvmix_zero)) +                &
            sqrt(max(Nsqr_iface(kt+1),cvmix_zero))) *             &
            0.5_cvmix_r8

Also, we want the default value for lavg_N_or_Nsqr to be .true. rather than .false.

The cvmix_put_kpp_real routine has a different structure from every other cvmix_put_[module] routine... it uses an optional cvmix_kpp_params_type argument and otherwise allows the user to edit variables in the private cvmix_kpp_params_saved variable. We want these routines to be consistent from module to module, and I like the optional argument structure.

Each module should have a private saved cvmix_[module]_params_type, and each public routine that takes that type as an argument should do so optionally - if an argument is present, it provides whatever parameters are needed otherwise the module-wide structure is used.

This option currently exists in

• cvmix_background
• cvmix_convection
• cvmix_ddiff
• cvmix_kpp

but is missing in

• cvmix_shear
• cvmix_tidal

We've talked via email about adding an error code as an optional argument (to allow modules to handle errors on their own rather than getting a stop call which is sure to cause issues with MPI). There would also be a cvmix_error module in src/ to handle the errors for the stand-alone drivers... likely just printing a statement followed by a stop call.

On the same note, we should put some thought into handling warnings - I think we should print warnings in a variety of cases (over-writing existing parameters, using cvmix_put to write a scalar into an array, etc)... but there should be an easy way to suppress the warnings. I'm leaning towards a namelist option that sets the global variable lwarn, but I'm open to other ideas - a namelist variable stored somewhere else? Copious ifdefs along with a compile-time setting? Going with the namelist option, I'm picturing a cvmix_warn routine that takes a string as an argument and prints "WARNING: " // warn_string... with the check against lwarn or whatever in the cvmix_warn routine. Maybe something like

module cvmix_warn

public :: cvmix_warn
public :: cvmix_set_lwarn

logical, private, save :: lwarn = .true.

...

end module cvmix_warn

So by default we print all warnings, but the user can call cvmix_set_lwarn(.false.) to suppress them (and then later call cvmix_set_lwarn(.true.) to turn them back on).

Right now, POP+CVMix can only run with a single thread per MPI task (by default, the NCAR machine runs POP with 2 threads per MPI task). At this point, I don't know if the issue is on the POP side or the CVMix side.

Low-level routines (cvmix_[module]_coeffs) should have separate inputs and outputs for the diffusivities so that the user can determine what to do with new diffusivities (overwrite old, add to old, etc). The interfaces using the CVMix wrapper should have some parameter flag to determine what to return (see cvmix_background and the 'handle_old_vals' parameter).

When computing the bulk Richardson number, we want

Rib = (d-dr)*(Br-B(d))/( |Vr-V(d)|^2 + Vt^2(d) )

Where d is the current cell-center depth and dr is the center of the surface layer. As @toddringler pointed out, dr = 0.5_surf_layer_ext_d, so instead of adding an additional argument to cvmix_kpp_compute_bulk_Richardson() we'll compute

Rib = (1-0.5_surf_layer_ext)_d*(Br-B(d))/( |Vr-V(d)|^2 + Vt^2(d) )

Need physical parameters for testing KPP. Ideally, would have multiple sets of surface buoyancy force, surface friction velocity, OBL_depths, and interior mixing values to exercise as many branches of the KPP code as possible. Specifically, for computing turbulent velocity scales, would like to test the case where surf_fric_vel = 0 & surf_buoy_force < 0 as well other general cases.

The Pacanowski-Philander (PP) version of shear mixing is based on the gradient Richardson number (RICH). The scheme proposed in PP (1981) clearly assumes RICH >= 0. But within CVMix, RICH is not tested to be 0 or positive when using PP. Note that when using shear KPP, RICH<0 is tested, so CVMix is assuming RICH can be either positive or negative.

While updating cvmix_shear to separate inputs and outputs / have a low-level and a wrapped version of cvmix_coeffs_shear, I don't know that I handled the "no_diff" operation correctly. I'll go back and check once the rest of the modules have separate inputs / outputs and low-level coeff routines.

@StephenGriffies pointed out two parameters that might be useful:

1. A minimum OBL_depth

2. A prescribed OBL_depth -- this could be helpful in diagnosing differences in the WSwPSBF tests

Per a discussion on September 23, 2013 -- the cvmix_put and cvmix_get routines improve the readability of the code but result in a performance hit. Is that a worthwhile trade-off?

Need to make sure parameters are set with reasonable values (in cvmix_[module]_init routines?). This is already done in some situations: kpp makes sure c_s and c_m are positive and that surf_layer_ext is in [0,1].

We cleaned up KPP before the beta release, but other modules still have allocatable arrays that should have a declared size...

After first public release, we should look into vectorizing the code. (This is a pretty vague ticket, details will be added via comments.)

See https://github.com/CVMix/CVMix-src/tree/feature/kpp_optimized_for_openmp_no_allocates for tips on removing allocate calls that slow down the code (GFDL reported 30% speed up using this branch rather than corresponding trunk)

Many cvmix_kpp routines were made public for testing, but really should only be callable from the cvmix_coeffs_kpp() routine [shape functions, kOBL_depth, enhanced_diff, compute_nonlocal, compute_nu_at_OBL_depth]

Instead of computing the enhanced diffusivity when using "MatchBoth" (and no other time), introduce a new parameter [cvmix_kpp_params%enhance_diff = .true. / .false.]. For now, we can default to enhance_diff = .true. but that may be changed later (or depending on choice of MatchTechnique).

If you want to write a variable to netCDF, you should be able to specify one of many names (for example, "SqrBuoyancy", "SqrBuoyancyFreq", "SqrBuoyancyFreq_iface", "buoyancy_iface", etc) but those should all result in a variable with the same name ("SqrBuoyancyFreq"?) being made in the netCDF file. Similarly, you should be able to edit attributes using a wide range of names.

I think the best way to do this will be to create a simplify_name routine so that the multiple names only need to be handled in a single routine.

Shear mixing will error out if you do not set all the necessary parameters while every other module sets the parameters to default values if they are not provided by the user.

Nothing in the reg_test/ directory exercises code in cvmix_convection module.

We should make sure the shape of any arrays being stored via cvmix_put are the same shape as the data type attribute they are being stored in.

Example:

1. This is good, because zt should be size nlev

allocate(depth(5))
cvmix_put(CVMix_vars, 'nlev', 5)
cvmix_put(CVMix_vars, 'zt', depth)

2. This is bad, but results in a seg-fault instead of a clean error

allocate(depth(8))
cvmix_put(CVMix_vars, 'nlev', 5)
cvmix_put(CVMix_vars, 'zt', depth)

3. This is also bad, and I think will also seg-fault... at the very least it will do something unexpected

allocate(depth(5))
cvmix_put(CVMix_vars, 'nlev', 8)
cvmix_put(CVMix_vars, 'zt', depth)

Need a clean way to have CVMix optionally output diagnostic variables.

Need a general cvmix_finalize to deallocate memory (or clean up pointers?) as well as method-specific cvmix_MODULE_finalize routines for methods that allocate memory (background mixing is the only one for now).

The variable naming convention should be the following:

1. Variables consisting of multiple words should capitalize the first letter of each word
2. Variables that are located at cell centers should have _cntr appended
3. Variables that are located at cell interfaces should have _iface appended
4. Declaration of variable should include full name and units (note that CVMix uses the 'mks' unit convention)

Exceptions to the above - diffusivities will be slightly abbreviated:

• Udiff_iface instead of MomentumDiffusivity_iface
• Tdiff_iface instead of TemperatureDiffusivity_iface
• Sdiff_iface instead of SalinityAndOtherTracerDiffusivity_iface

Same standards should be applied to intent(in)s and intent(out)s.

The Simmons scheme uses a static tidal energy field, and the vertical deposition is also static in time, so we only need to compute energy_flux*vert_dep/rho_fw once (rather than every time step, as is currently the case). I'd like to refactor the tidal mixing module to include a compute_tidal_diffusivity() routine that can be during the initialization phase of the GCM to compute this value and store it in the cvmix_data_type structure, then cvmix_coeffs_tidal_low() can reference the stored data rather than computing it every time.

Testing should be split into two directories:

1. regression tests -- basically all of the tests we have now are regression tests
2. unit tests -- testing at a very basic level. Do things like ensure cvmix_put and cvmix_get write to the proper place, the initialization routines provide the expected default values, etc etc.

Also, should provide some light-weight examples for both stand-alone single column cases and (if possible) an outline of steps needed to import into a pre-existing ocean model.

Random thoughts:

1. Examples could build libcvmix.a out of src/shared/ and then link a simple driver from each subdirectory
2. src/drivers/ could be thought of as src/reg_tests/
3. What's the best way to build unit tests? All in single Fortran file (that, like examples, would only depend on libcvmix.a)? Should we require pfunit?
4. We should really look into automating the regression tests across a handful of machines / compilers.

static_diff(:,2) and static_visc(:) don't mesh with the renamed Tdiff, Sdiff, and Mdiff.

Shear mixing needs two background mixing parameters for Pacanowski-Philander. There are three options for how to handle this:

1. Add background parameters to the shear mixing parameter type
2. Add a variable of background parameter type to the shear mixing parameter type
3. Have the user set a background parameter variable type and pass it to the shear module when using PP

Currently option (3) is in use, but I think it's the worst option of the group. I'm not sure which of the first two would be the best to use, though. I'm leaning towards (2) because that already has all the needed support for 1D or 2D variations in the background field.

The way CVMix handles the nonlocal term doesn't account for enhanced diffusivity... need to multiply by Kenh/Kold at bottom of boundary layer

There have been a lot of changes to cvmix_coeffs_kpp_low between June 18 (4d2d987) and June 27 (61292f1), and it has changed answers in the GFDL single-column test case... Before making the beta tag, we want to make sure we understand what is happening in the routine with the GFDL default options and that everything looks okay.

Currently cvmix_kpp is the only module that allows users to pass in arrays rather than requiring users to pack data in to the CVmix_data_type structure. The current cvmix_[module]coeffs should be renamed cvmix[module]coeffs_wrap and it should call cvmix[module]_coeffs_low, where computation is actually done.

v0.53:

cvmix_kinds_and_types.F90 says:
! For KPP we need buoyancy (as opposed to buoyancy frequency)

then defines:
real(cvmix_r8), dimension(:), pointer :: buoyancy_cntr => NULL() ! units: m s^-2

but I can not find where/how buoyancy_cntr is used in the code. Is buoyancy_cntr used in CVMix/KPP?

a "grep buoyancy_cntr *.F90" turns up nothing except puts/gets.

As I have been looking through CVMIX, I have an issue with the following code in cvmix_kpp_compute_turbulent_scales_1d. In it, zeta (sigma * h / L) is computed as follows

 if (surf_fric_vel.ne.cvmix_zero) then
      do kw=1,n_sigma
        ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
        ! at surf_layer_ext
        zeta(kw) = min(surf_layer_ext, sigma_coord(kw)) * OBL_depth *         &
                   surf_buoy_force*vonkar/(surf_fric_vel**3)
      end do

I am wondering why the convective case limitation on sigma is included for stable buoyancy forcing as well? Large et al. (1994), Troen and Mahrt (1986), and the CVMIX description document assume that the similarity functions for the surface layer extend from 0 < sigma < 1

I have run a test case with positive buoyancy forcing and a constant wind stress in a SCM framework. Below is a plot of temperature profiles through time, with the kpp diagnosed boundary layer depth in black.

The corresponding result from a similarly forced LES is below (the line is a diagnosed mixed layer depth from a temperature threshold criterion)

I have then run a test that changed the CVMIX code above to

    if (surf_fric_vel.ne.cvmix_zero) then
      do kw=1,n_sigma
        ! compute scales at sigma if sigma < surf_layer_ext, otherwise compute
        ! at surf_layer_ext
        if(surf_buoy_force .ge. cvmix_zero) then
                zeta(kw) = sigma_coord(kw) * OBL_depth *                      &
                        surf_buoy_force*vonkar/(surf_fric_vel**3)
        else
                zeta(kw) = min(surf_layer_ext, sigma_coord(kw)) * OBL_depth *         &
                   surf_buoy_force*vonkar/(surf_fric_vel**3)
        endif
      end do

with this change, the temperature profile plot becomes

which is closer to the LES result.

It seems that the in stable conditions, sigma should not be confined less than the surface layer extent.

Obtained with intel, g77 and pgi compilers using -g. I suspect that OBL_Mdiff, OBL_Tdiff and OBL_Sdiff should be allocated with one extra slot.

Image              PC                Routine            Line        Source             
MOM6               00000000097E558A  cvmix_kpp_mp_cvmi         797  cvmix_kpp.F90
MOM6               000000000689F111  mom_kpp_mp_kpp_ca         617  MOM_KPP.F90
MOM6               0000000008D6016C  mom_diabatic_driv         473  MOM_diabatic_driver.F90
MOM6               000000000B99E375  mom_mp_step_mom_          803  MOM.F90
MOM6               0000000009B52933  ocean_model_mod_m         368  ocean_model_MOM.F90
MOM6               000000000154C038  MAIN__                    695  coupler_main.F90
forrtl: severe (408): fort: (2): Subscript #1 of the array OBL_MDIFF has value 6 which is greater than the upper bound of 5
'''

The non-local transport returned from cvmix_coeffs_kpp_low() contains a NaN in the first element (surface). Bug appeared in commit 4d2d987. Bug is still manifest in 4f3943a.

A possibly related change is that the diffusivity at the surface is no longer zero but a smallish (not tiny 1e-5) value.

Configuration uses MatchTechnique='SimpleShapes'.

The last working commit was 7ec9961.

I think the code is more readable if all routines are written in the same style, and the put routines exist for a reason.

Per a discussion on September 23, 2013 -- should we add new interface for compute_turbulent_shear() to allow for array of OBL_depth (used for computing Ri_bulk)? In more general terms, should more routines have interfaces to handle multiple columns at once?

POP's convective mixing scheme is only applied below the boundary layer. I have a branch of the repository (on my personal fork) where I've added a new option to the convective mixing initialization named lnoOBL, with a default value of .true. implying we do not apply in the boundary layer. I don't know if

(a) this is a feature we want in CVMix (alternatively, GCMs can zero out the convective mixing contribution in the OBL), or
(b) if this is the best implementation.

Some shortfallings -- currently, lnoOBL is only checked if we are using Brunt-Vaisala convective mixing with BVsqr_convect >= 0. Once I have a better idea on the desired application / defaults / variable names I'll clean up the code and apply it regardless of convective mixing scheme.

Add parameter for min_Vtsqr (default 1e-10 m^2/s^2), have cvmix_kpp_compute_unresolved_shear return

min(computed Vt^2, min_Vtsqr)

In LMD94, the nonlocal term is defined (by manipulating Eqs (10) and (20) by

K*gamma = -Cs * G(sigma) * Q

(the minus sign is the CVMix convention and CVMix returns K*gamma/Q to the base model).

@StephenGriffies pointed out that LMD assumes G(0) = 0 => the nonlocal term is 0 at the surface... but it may be desirable to use a shape function that allows K*gamma = -Q at the surface.

I think there are two ways to accomplish this. Both start with removing cvmix_kpp_compute_nonlocal() [which returns G(sigma)*Cs] and simply computing G(sigma). From there, we can either:

1. As Steve and I discussed over the phone, we can add Cs to the cvmix_kpp_params_type. When you call cvmix_init_kpp(), if you omit the Cs argument then it computes Cs = C* * kappa * (__)^1/3 or you can set it to 1

2. Something else that I just thought of would be to create Cs as a module-wide parameter. If you choose MatchTechnique = 'ParabolicNonlocal' then Cs = 1, otherwise Cs = the value above.

In either case, you multiply Tnonlocal and Snonlocal by Cs after computing all the G(sigma) values in the column.

Relatedly, Steve has suggested two other shape functions for the nonlocal term, a linear function and a cubic function. I think the best way to implement those would be to have another MatchTechnique parameter (so MatchTechnique controls the local shape function and MatchTechnique2 or MatchTechniqueNonlocal controls the nonlocal shape function). I'll tackle that after fixing this issue with Cs, regardless of whether we go with option (1) or option (2) it should be a straight-forward update.

The original PP scheme includes a background term that is added on.

The implementation in CVMix is the same and, as a result, breaks the "a al carte" approach that we strive for in CVMix. In order to run the CVMix shear PP module we have to initialize CVMix background module, i.e. shear depends on background.

I propose that we drop the background part of PP. Most everyone who uses CVMix shear will be using some form of background mixing that is set outside of PP. For example, right now MPAS-O has to test if shear mixing is using PP and, if so, subtract out the background from the result (because we have added it elsewhere.)

It will lead to clearer code all around to simply remove the addition of background from CVMix shear PP. (Note that it is not added into CVMix shear KPP).

In cvmix_kpp_coeffs_low, the indexing in variables like Gat1, DGat1, visc_at_OBL isn't very clear... I think I might put in a module-wide parameter index (CVMIX_KPP_MDIFF_IND, etc) to clarify.

We need to use N / N^2 at the bottom interface of a cell when computing unresolved shear because the parameter Cv (and the term Cv*N) is defined in this manner. Computing N at the cell center would require recomputing Cv.

I believe cvmix_kpp requires Sdiff at this time, but the actual KPP routine doesn't require double diffusion and therefore flops / memory can be saved for users who only have momentum and tracer diffusion coefficients.

A few functions within the CVMix KPP interface take the argument max_nlev to allow KPP to be computed over a subset of the entire column (e.g. cvmix_coeffs_kpp).

1. This option should be available to all routines (e.g. compute_turbulent_scales, compute_bulk_richardson).

2. Even when the max_nlev option is included in a routine, there are portions of that routine that still compute over nlevs For example, within coeffs_kpp, compute_turbulent_scales is called it is over nlevs (determined by the size of zw_ctr). Thus, if max_nlev != nlevs and all of zw_ctr is sent, an issue can arise.

The Makefiles I've written get the job done, though occasionally lack of maintenance has led to issues with some of the compilers / targets. It would probably be cleaner and easier to manage if we switched to something like CMake

Need clear instructions for members of community that want to add mixing techniques to the code. For example:

1. Design requirements for cvmix_[module].F90
2. Stand-alone driver / regression test
3. Documentation (both in-code and a chapter for the manual)

We want each module to have a clear owner so that users can ask the author directly instead of going through us.

I only ran into this issue when using KPP, but the problem ends up being a combination of lines:
451 and 600-602 in https://github.com/CVMix/CVMix-src/blob/master/src/shared/cvmix_kpp.F90

Basically, the arrays new_Mdiff (and such) are allocated using nlev, which is the number of active levels in a column, not necessarily the number of total levels in the column. Later a direct copy happens moving old_Mdiff into Mdiff_out, but they have different sizes (one is active, one is total) which causes an array out of bounds with some compilers. Ideally, the arrays would be allocated using the total number of levels whether or not they are used.

Obviously some system is needed to ensure code on the master branch compiles and runs as expected... ideally with little or no human intervention.

Necessary Features:

1. Build executable with all compilers available on machine, report build errors
2. Run regression tests from reg_tests/ directory, report run errors
3. Compare output to baseline output, report answer changes

Some discussion has taken place on the developer's mailing list:

• Recommendation to avoid shell scripts as much as possible
• Look into CMake / CTest?

Besides needing to create this testing infrastructure, I think we need to completely overhaul reg_tests/ so that it will play nice with the test system and be more amenable to introducing new regression tests.