Last edited Timestamp?

**Author** : Matthieu Leclair

**ticket** : #863

**Branch** : z_tilde coordinate

**Secondary author** : Andrew Coward

**Branch** : reimplementation of z_tilde coordinate starting from 3.5alpha

wiki:ticket/863/reimplementation?

## The z_tilde vertical coordinate

Implementation of the z_tilde Arbitrary Lagrangian Eulerian vertical coordinate discribed in Leclair & Madec, Ocean Modelling, 2011

Baroclinic high frequency vertical motions are followed by the coordinate (Lagrangian part) while barotropic and baroclinic low frequency ones are discribed through cross-level velocities (Eulerian part).

## Structure

z_star and z_tilde vertical coordinates are now two possible "sub-options" of the vvl key. A third one, the full layer coordinate has been introduced for academic and debuging purposes.

Next sea surface height is first computed. The latter is necessary before computing next vertical scale factors. And now vertical velocity (eulerian cross-level velocity) can then be calculated.

After this sequence, tracers and dynamics terms can be computed. Next tracers and dynamics fields are updated.

Finally, after the swap of those fields, sea surface height and vertical scale factors can be swaped.

Calling sequence in `step` :

!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Update data, open boundaries, surface boundary condition (including sea-ice) !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< . . . !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Ocean dynamics : hdiv, rot, ssh, e3, wn !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CALL ssh_nxt ( kstp ) ! after ssh IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors CALL wzv ( kstp ) ! now cross-level velocity !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Ocean physics update (ua, va, ta, sa used as workspace) !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< . . . !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Passive Tracer Model !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< . . . !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Active tracers (ua, va used as workspace) !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< . . . !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Dynamics (ta, sa used as workspace) !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< . . . CALL ssh_swp( kstp ) ! swap of sea surface height IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors

## Newly created , renamed and suppressed variables

- 2D fields
`sshu_b`,`sshu_n`,`sshu_a`,`sshv_b`,`sshv_n`,`sshv_a`and`sshf_n`needed for the previous z_star implementation are suppressed. - 3D fields
`mut`,`muu`,`muv`and`muf`and 2D fields`ee_t`,`ee_u`,`ee_v`,`ee_w`needed for the previous z_star implementation are also suppressed. `e3t_1`,`e3u_1`,`e3v_1`,`e3w_1`,`e3f_1`,`gdep3w_1`,`gdept_1`and`gdepw_1`are renamed with`_n`instead of`_1.`- 3D fields
`e3t_b`,`e3u_b`,`e3v_b`,`e3uw_b`,`e3vw_b`and`e3t_a`,`e3u_a`,`e3v_a`are created. - 3D fields
`e3t_t_b`,`e3t_t_n`and`e3t_t_a`(baroclinic part of the scale factors in z_tilde case) are created. - 3D field
`hdiv_lf`(low frequency part of the baroclinic horizontal divergence) is created. - 2D fields
`frq_restore_e3t`and f`rq_restore_hdv (restoring frequencies) are created.` - 2D fields
`e12t`,`e12t_1`,`e12u`,`e12u_1`,`e12v`,`e12v_1`,`e12f`,`e12f_1`have been created.

All of them are not used but it will be usefull to use them in the future in many routines

Budget: 8 3D fields added and 1 2D field removed

renamed :

- 1D fields
`gdept_0`,`gdepw_0,...`,`e3t_0`,`e3w_0`, …. etc have been renamed with suffix _1d instead of _0. (see Changeset 3839)

## Core description

### Step 1 : Initialisation

`dom_vvl_init` is called by `domain:`

- The namelist nam_vvl is read and options are controlled. there are 3 possible variable vertical coordinates

- z_star coordinate: spread the ssh increment over the water column at each time step.

It is equivalent to a Lagrangian treatment of the barotropic flow. - z_tilde coordinate: Treat barotropic flow in the same way as z_star and also absorb high frequency baroclinic motions in the vertical scale factor variations.
- layer: absorb all vertical motions with the variations of vertical scale factors.

- Allocate arrays depending on the chosen options.

- Read restart file if needed.

- Initialize before and now scale factors at different grid points.

Scale factor anomalies at T-points are interpolated to another grid points and added to the background scale factor at this point. The interpolation is a surface weighted averaging in the horizontal direction and simple half averaging in the vertical direction

### Step 2 : Next sea surface height

Next sea surface height is computed by `ssh_nxt` routine in module `sshwzv`.

### Step 3 : Next vertical scale factor

Next vertical scale factors are computed by `dom_vvl_nxt` routine in `domvvl` module. 2 cases are distinguished:

**z_star coordinate**

Vertical scale factors are not obtained with the repartition of the sea surface heigt over the vertical levels anymore. The reason is that the conservation correction applied to the leap-fog time stepping scheme for tracers and volume only concerns the first ocean level. We thus have to treat the entire 3D vertical scale factors:

z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) DO jk = 1, jpkm1 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) END DO

**z_tilde or layer coordinate**

1. Baroclinic Vertical scale factor anomalies (`e3t_t_a`) are obtained thanks to 3 tendency terms:

* High frequency baroclinic horizontal divergence (which is the entire divergence`hdivn`in the layer case).

* Restoring towards 0 (only for z_tilde).

* Thickness diffusion term.

2. Barotropic part of the scale factor is incremented by the repartition of the sea surface height variation (as in the z_star case)

### Step 4 : Now vertical velocity (Eulerian cross-level velocity)

Now Eulerian vertical velocity (`wn`) is computed by `wzv` routine in `sshwzv` module. There is no particular reason for letting this routine in this module but no real reason to put it elsewhere … The computation of wn takes the thickness diffusion transports (`un_td` and `vn_td`) into account in order to avoid a compensating vertical velocity.

### Step 5 : Thickness diffusion velocity in tracers advection

In order to be consistent between volume and tracer equations the thickness diffusion transport is also taken into account in the lateral tracer advection terms.

### Step 6 : Swap of sea surface height

This part was previously done in `ssh_nxt` routine which could cause a bit of confusion here. the routine is now called ssh_swp, still in the `sshwzv` module.

### Step 7 : Swap of vertical scale factors

Swap of vertical scale factors is done in two different routines.

- The first part concerns
`fse3t_b`,`fse3u_b`and`fse3v_b`. This is done in`dynnxt`routine because they are needed there. - The second part concerns
`fse3t_n`,`fse3u_n`and`fse3v_n`. This is done by`dom_vvl_swp`routine in`domvvl`module. - The third part concerns
`e3t_t_b`and`e3t_t_n`. This is done by`dom_vvl_swp`routine in`domvvl`module.

It could be done in`dom_vvl_sf_nxt`which would make`e3t_t_a`purely local to this routine. A simple work array could thus be used (e.g. ze3t_t_a) instead of a saved module array. In this case the entire module would be less readable

`dom_vvl_swp` routine also interpolates scale factors from t-, u- or v-points to other grid points and computes depths and water column heights.

Write restart and outputs.

## Side modifications

### Time splitting external mode

Modifications in the vvl case are done. The previous version only treated the z_star option. It is now more general.

### Coupling between ocean and sea ice

In the vvl case `fse3t_m` is now a statement function for `e3t_m` whis is a 2D variable. In non-vvl cases `fse3t_m(:,:)` stands for fse3t_0(:,:,1).`fse3t_m` is treated in the vvl case as all other `*_m` variables in `sbcssm` module.

### Tracer lateral diffusion

The 2D fields `e2_1u = e2u / e1u` and `e1_2v = e1v / e2v` are used in the computation of the thickness diffusion term. They were allready computed and saved locally in `traldf_lap` module under the names e1ur and e2vr. Because they are used in all other tracer diffusion routines (and now the thickness diffusion term in `dom_vvl_sf_nxt`), they are now declared in `dom_oce` and computed in `domain` just after all e12x and e12x_1 fields. The different tracer lateral diffusion routines have thus been modified to use these new variables. This could be done in many other routines.

### Isopycnal diffusion

The computation of potential density is done at a certain depth. In the vvl case, neighbor cells potential densities are not computed at the exact same depth, which causes an error in the computation of isopycnal slopes. To fix this problem, the depth has been added as an argument in the `eos_insitu` and `eos_insitu_pot` routines.

When the `eos` interface is called before computing the isopycnal slopes, `fsdept_0` is used. When it's called before the computation of the hydrostatic pressure gradient, the current depth `fsdept_n` is used.

### Thickness weighted outputs (active and passive tracers + velocities)

In vvl case, the "toce" variable in the netcdf file does not contain temperature anymore but the heat content fse3t_n * tn.

The same thing is done for salinity and horizontal velocities (only in flux form advection but I'm not sure if it has to be done in any vvl case).

The same treatment is also applied to passive tracers.

3 other variables are created in netcdf outputs:

`e3t_n`is the vertical scale factor at t-points (`_n`suffix indicates that it is different from the background scale factor).

e3u and e3v and e3w are not stored, they can be recomputed off-line.`dept_n`is the depth at t-points with a 0 reference when the ocean is at rest (ssh is removed).

Depths at u, v and w points should be computed off-line also.`e3tdef`is the squared deformation ratio of`e3t`given in %^2.

During the post treatment, the square root of this quantity gives the variance of the level thickness in %.

z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - fse3t_0(:,:,:) ) / fse3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 CALL iom_put( "e3t_n" , fse3t_n (:,:,:) ) CALL iom_put( "dept_n" , fsde3w_n (:,:,:) ) CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )

## To be done

### Compatibility with een vorticity scheme

The vertical scale factor at f points should remain constant in time if the f-point belongs to an "earth cell".

IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t) DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, jpim1 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25 IF( ze3f(ji,jj,jk) /= 0._wp ) ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk) END DO END DO END DO CALL lbc_lnk( ze3f, 'F', 1. ) ENDIF

## Testing

Testing could consider (where appropriate) other configurations in addition to NVTK].

NVTK Tested | '''YES/NO''' |

Other model configurations | '''YES/NO''' |

Processor configurations tested | [ Enter processor configs tested here ] |

If adding new functionality please confirm that the New code doesn't change results when it is switched off and ''works'' when switched on | '''YES/NO/NA''' |

(Answering UNSURE is likely to generate further questions from reviewers.)

'Please add further summary details here'

- Processor configurations tested
- etc——

## Bit Comparability

Does this change preserve answers in your tested standard configurations (to the last bit) ? | '''YES/NO ''' |

Does this change bit compare across various processor configurations. (1xM, Nx1 and MxN are recommended) | '''YES/NO''' |

Is this change expected to preserve answers in all possible model configurations? | '''YES/NO''' |

Is this change expected to preserve all diagnostics? ,,''Preserving answers in model runs does not necessarily imply preserved diagnostics. '' | '''YES/NO''' |

If you answered '''NO''' to any of the above, please provide further details:

- Which routine(s) are causing the difference?
- Why the changes are not protected by a logical switch or new section-version
- What is needed to achieve regression with the previous model release (e.g. a regression branch, hand-edits etc). If this is not possible, explain why not.
- What do you expect to see occur in the test harness jobs?
- Which diagnostics have you altered and why have they changed?Please add details here……..

## System Changes

Does your change alter namelists? | '''YES/NO ''' |

Does your change require a change in compiler options? | '''YES/NO ''' |

If any of these apply, please document the changes required here…….

## Resources

''Please ''summarize'' any changes in runtime or memory use caused by this change……''

## IPR issues

Has the code been wholly (100%) produced by NEMO developers staff working exclusively on NEMO? | '''YES/ NO ''' |

If No:

- Identify the collaboration agreement details
- Ensure the code routine header is in accordance with the agreement, (Copyright/Redistribution? etc).Add further details here if required……….