| Version 26 (modified by mlelod, 12 years ago) (diff) |
|---|
Last edited Timestamp?
Author : Matthieu Leclair
ticket : #863
Branch : z_tilde coordinate
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 frq_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
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 and fse3t_n and is done in dynnxt routine because they are needed there.
- The second part concerns e3t_t_b and e3t_t_n. This is done by dom_vvl_swp routine in domvvl module.
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(:,:,:) )
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..........
