Context Navigation

domvvl.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 55.6 KB

Line
1	MODULE domvvl
2	!!======================================================================
3	!! * MODULE domvvl *
4	!! Ocean :
5	!!======================================================================
6	!! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code
7	!! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate
8	!! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
9	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
10	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
14	!! dom_vvl_init : define initial vertical scale factors, depths and column thickness
15	!! dom_vvl_sf_nxt : Compute next vertical scale factors
16	!! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid
17	!! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
18	!! dom_vvl_rst : read/write restart file
19	!! dom_vvl_ctl : Check the vvl options
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers
22	USE phycst ! physical constant
23	USE dom_oce ! ocean space and time domain
24	USE sbc_oce ! ocean surface boundary condition
25	USE wet_dry ! wetting and drying
26	USE usrdef_istate ! user defined initial state (wad only)
27	USE restart ! ocean restart
28	!
29	USE in_out_manager ! I/O manager
30	USE iom ! I/O manager library
31	USE lib_mpp ! distributed memory computing library
32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
33	USE timing ! Timing
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC dom_vvl_init ! called by domain.F90
39	PUBLIC dom_vvl_zgr ! called by isfcpl.F90
40	PUBLIC dom_vvl_sf_nxt ! called by step.F90
41	PUBLIC dom_vvl_sf_update ! called by step.F90
42	PUBLIC dom_vvl_interpol ! called by dynnxt.F90
43
44	! !!* Namelist nam_vvl
45	LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
46	LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate
47	LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate
48	LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate
49	LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate
50	LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer
51	! ! conservation: not used yet
52	REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient
53	REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days]
54	REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days]
55	REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation
56	LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints
57
58	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport
59	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence
60	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors
61	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors
62	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors
63	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence
64
65	!! * Substitutions
66	# include "do_loop_substitute.h90"
67	!!----------------------------------------------------------------------
68	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
69	!! $Id$
70	!! Software governed by the CeCILL license (see ./LICENSE)
71	!!----------------------------------------------------------------------
72	CONTAINS
73
74	INTEGER FUNCTION dom_vvl_alloc()
75	!!----------------------------------------------------------------------
76	!! * FUNCTION dom_vvl_alloc *
77	!!----------------------------------------------------------------------
78	IF( ln_vvl_zstar ) dom_vvl_alloc = 0
79	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
80	ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , &
81	& dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , &
82	& STAT = dom_vvl_alloc )
83	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
84	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
85	un_td = 0._wp
86	vn_td = 0._wp
87	ENDIF
88	IF( ln_vvl_ztilde ) THEN
89	ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
90	CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
91	IF( dom_vvl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
92	ENDIF
93	!
94	END FUNCTION dom_vvl_alloc
95
96
97	SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
98	!!----------------------------------------------------------------------
99	!! * ROUTINE dom_vvl_init *
100	!!
101	!! ** Purpose : Initialization of all scale factors, depths
102	!! and water column heights
103	!!
104	!! ** Method : - use restart file and/or initialize
105	!! - interpolate scale factors
106	!!
107	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
108	!! - Regrid: e3[u/v](:,:,:,Kmm)
109	!! e3[u/v](:,:,:,Kmm)
110	!! e3w(:,:,:,Kmm)
111	!! e3[u/v]w_b
112	!! e3[u/v]w_n
113	!! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
114	!! - h(t/u/v)_0
115	!! - frq_rst_e3t and frq_rst_hdv
116	!!
117	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
118	!!----------------------------------------------------------------------
119	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
120	!
121	IF(lwp) WRITE(numout,*)
122	IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
123	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
124	!
125	CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)
126	!
127	! ! Allocate module arrays
128	IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
129	!
130	! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
131	CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
132	e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all
133	!
134	CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
135	!
136	END SUBROUTINE dom_vvl_init
137	!
138	SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa)
139	!!----------------------------------------------------------------------
140	!! * ROUTINE dom_vvl_init *
141	!!
142	!! ** Purpose : Interpolation of all scale factors,
143	!! depths and water column heights
144	!!
145	!! ** Method : - interpolate scale factors
146	!!
147	!! ** Action : - e3t_(n/b) and tilde_e3t_(n/b)
148	!! - Regrid: e3(u/v)_n
149	!! e3(u/v)_b
150	!! e3w_n
151	!! e3(u/v)w_b
152	!! e3(u/v)w_n
153	!! gdept_n, gdepw_n and gde3w_n
154	!! - h(t/u/v)_0
155	!! - frq_rst_e3t and frq_rst_hdv
156	!!
157	!! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling.
158	!!----------------------------------------------------------------------
159	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
160	!!----------------------------------------------------------------------
161	INTEGER :: ji, jj, jk
162	INTEGER :: ii0, ii1, ij0, ij1
163	REAL(wp):: zcoef
164	!!----------------------------------------------------------------------
165	!
166	! !== Set of all other vertical scale factors ==! (now and before)
167	! ! Horizontal interpolation of e3t
168	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U
169	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
170	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V
171	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
172	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F
173	! ! Vertical interpolation of e3t,u,v
174	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W
175	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' )
176	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW
177	CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
178	CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW
179	CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
180
181	! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
182	e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
183	e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
184	e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
185	!
186	! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep)
187	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration)
188	gdepw(:,:,1,Kmm) = 0.0_wp
189	gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg
190	gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
191	gdepw(:,:,1,Kbb) = 0.0_wp
192	DO_3D_11_11( 2, jpk )
193	! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
194	! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
195	! ! 0.5 where jk = mikt
196	!!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ??
197	zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
198	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
199	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) &
200	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm))
201	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
202	gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
203	gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) &
204	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb))
205	END_3D
206	!
207	! !== thickness of the water column !! (ocean portion only)
208	ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....
209	hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
210	hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
211	hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
212	hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
213	DO jk = 2, jpkm1
214	ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
215	hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
216	hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
217	hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
218	hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
219	END DO
220	!
221	! !== inverse of water column thickness ==! (u- and v- points)
222	r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF
223	r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
224	r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
225	r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
226
227	! !== z_tilde coordinate case ==! (Restoring frequencies)
228	IF( ln_vvl_ztilde ) THEN
229	!!gm : idea: add here a READ in a file of custumized restoring frequency
230	! ! Values in days provided via the namelist
231	! ! use rsmall to avoid possible division by zero errors with faulty settings
232	frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp )
233	frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
234	!
235	IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile
236	frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings
237	frq_rst_hdv(:,:) = 1._wp / rdt
238	ENDIF
239	IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator
240	DO_2D_11_11
241	!!gm case \|gphi\| >= 6 degrees is useless initialized just above by default
242	IF( ABS(gphit(ji,jj)) >= 6.) THEN
243	! values outside the equatorial band and transition zone (ztilde)
244	frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp )
245	frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
246	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star
247	! values inside the equatorial band (ztilde as zstar)
248	frq_rst_e3t(ji,jj) = 0.0_wp
249	frq_rst_hdv(ji,jj) = 1.0_wp / rdt
250	ELSE ! transition band (2.5 to 6 degrees N/S)
251	! ! (linearly transition from z-tilde to z-star)
252	frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp &
253	& * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
254	& * 180._wp / 3.5_wp ) )
255	frq_rst_hdv(ji,jj) = (1.0_wp / rdt) &
256	& + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp &
257	& * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) &
258	& * 180._wp / 3.5_wp ) )
259	ENDIF
260	END_2D
261	IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
262	IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
263	ii0 = 103 ; ii1 = 111
264	ij0 = 128 ; ij1 = 135 ;
265	frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
266	frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt
267	ENDIF
268	ENDIF
269	ENDIF
270	ENDIF
271	!
272	IF(lwxios) THEN
273	! define variables in restart file when writing with XIOS
274	CALL iom_set_rstw_var_active('e3t_b')
275	CALL iom_set_rstw_var_active('e3t_n')
276	! ! ----------------------- !
277	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
278	! ! ----------------------- !
279	CALL iom_set_rstw_var_active('tilde_e3t_b')
280	CALL iom_set_rstw_var_active('tilde_e3t_n')
281	END IF
282	! ! -------------!
283	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
284	! ! ------------ !
285	CALL iom_set_rstw_var_active('hdiv_lf')
286	ENDIF
287	!
288	ENDIF
289	!
290	END SUBROUTINE dom_vvl_zgr
291
292
293	SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
294	!!----------------------------------------------------------------------
295	!! * ROUTINE dom_vvl_sf_nxt *
296	!!
297	!! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt,
298	!! tranxt and dynspg routines
299	!!
300	!! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness.
301	!! - z_tilde_case: after scale factor increment =
302	!! high frequency part of horizontal divergence
303	!! + retsoring towards the background grid
304	!! + thickness difusion
305	!! Then repartition of ssh INCREMENT proportionnaly
306	!! to the "baroclinic" level thickness.
307	!!
308	!! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case
309	!! - tilde_e3t_a: after increment of vertical scale factor
310	!! in z_tilde case
311	!! - e3(t/u/v)_a
312	!!
313	!! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling.
314	!!----------------------------------------------------------------------
315	INTEGER, INTENT( in ) :: kt ! time step
316	INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step
317	INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence
318	!
319	INTEGER :: ji, jj, jk ! dummy loop indices
320	INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers
321	REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars
322	LOGICAL :: ll_do_bclinic ! local logical
323	REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv
324	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t
325	!!----------------------------------------------------------------------
326	!
327	IF( ln_linssh ) RETURN ! No calculation in linear free surface
328	!
329	IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt')
330	!
331	IF( kt == nit000 ) THEN
332	IF(lwp) WRITE(numout,*)
333	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
334	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
335	ENDIF
336
337	ll_do_bclinic = .TRUE.
338	IF( PRESENT(kcall) ) THEN
339	IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
340	ENDIF
341
342	! ******************************* !
343	! After acale factors at t-points !
344	! ******************************* !
345	! ! --------------------------------------------- !
346	! ! z_star coordinate and barotropic z-tilde part !
347	! ! --------------------------------------------- !
348	!
349	z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
350	DO jk = 1, jpkm1
351	! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
352	e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
353	END DO
354	!
355	IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !
356	! ! ------baroclinic part------ !
357	! I - initialization
358	! ==================
359
360	! 1 - barotropic divergence
361	! -------------------------
362	zhdiv(:,:) = 0._wp
363	zht(:,:) = 0._wp
364	DO jk = 1, jpkm1
365	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
366	zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
367	END DO
368	zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
369
370	! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only)
371	! --------------------------------------------------
372	IF( ln_vvl_ztilde ) THEN
373	IF( kt > nit000 ) THEN
374	DO jk = 1, jpkm1
375	hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) &
376	& * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
377	END DO
378	ENDIF
379	ENDIF
380
381	! II - after z_tilde increments of vertical scale factors
382	! =======================================================
383	tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms
384
385	! 1 - High frequency divergence term
386	! ----------------------------------
387	IF( ln_vvl_ztilde ) THEN ! z_tilde case
388	DO jk = 1, jpkm1
389	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
390	END DO
391	ELSE ! layer case
392	DO jk = 1, jpkm1
393	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
394	END DO
395	ENDIF
396
397	! 2 - Restoring term (z-tilde case only)
398	! ------------------
399	IF( ln_vvl_ztilde ) THEN
400	DO jk = 1, jpk
401	tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
402	END DO
403	ENDIF
404
405	! 3 - Thickness diffusion term
406	! ----------------------------
407	zwu(:,:) = 0._wp
408	zwv(:,:) = 0._wp
409	DO_3D_10_10( 1, jpkm1 )
410	un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) &
411	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )
412	vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) &
413	& * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )
414	zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
415	zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
416	END_3D
417	DO_2D_11_11
418	un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
419	vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
420	END_2D
421	DO_3D_00_00( 1, jpkm1 )
422	tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &
423	& + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &
424	& ) * r1_e1e2t(ji,jj)
425	END_3D
426	! ! d - thickness diffusion transport: boundary conditions
427	! (stored for tracer advction and continuity equation)
428	CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
429
430	! 4 - Time stepping of baroclinic scale factors
431	! ---------------------------------------------
432	! Leapfrog time stepping
433	! ~~~~~~~~~~~~~~~~~~~~~~
434	IF( neuler == 0 .AND. kt == nit000 ) THEN
435	z2dt = rdt
436	ELSE
437	z2dt = 2.0_wp * rdt
438	ENDIF
439	CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
440	tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
441
442	! Maximum deformation control
443	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
444	ze3t(:,:,jpk) = 0._wp
445	DO jk = 1, jpkm1
446	ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
447	END DO
448	z_tmax = MAXVAL( ze3t(:,:,:) )
449	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
450	z_tmin = MINVAL( ze3t(:,:,:) )
451	CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain
452	! - ML - test: for the moment, stop simulation for too large e3_t variations
453	IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
454	IF( lk_mpp ) THEN
455	CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max )
456	CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min )
457	ELSE
458	ijk_max = MAXLOC( ze3t(:,:,:) )
459	ijk_max(1) = ijk_max(1) + nimpp - 1
460	ijk_max(2) = ijk_max(2) + njmpp - 1
461	ijk_min = MINLOC( ze3t(:,:,:) )
462	ijk_min(1) = ijk_min(1) + nimpp - 1
463	ijk_min(2) = ijk_min(2) + njmpp - 1
464	ENDIF
465	IF (lwp) THEN
466	WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
467	WRITE(numout, *) 'at i, j, k=', ijk_max
468	WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
469	WRITE(numout, *) 'at i, j, k=', ijk_min
470	CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
471	ENDIF
472	ENDIF
473	! - ML - end test
474	! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
475	tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )
476	tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
477
478	!
479	! "tilda" change in the after scale factor
480	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
481	DO jk = 1, jpkm1
482	dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
483	END DO
484	! III - Barotropic repartition of the sea surface height over the baroclinic profile
485	! ==================================================================================
486	! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
487	! - ML - baroclinicity error should be better treated in the future
488	! i.e. locally and not spread over the water column.
489	! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
490	zht(:,:) = 0.
491	DO jk = 1, jpkm1
492	zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
493	END DO
494	z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
495	DO jk = 1, jpkm1
496	dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
497	END DO
498
499	ENDIF
500
501	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !
502	! ! ---baroclinic part--------- !
503	DO jk = 1, jpkm1
504	e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
505	END DO
506	ENDIF
507
508	IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging
509	!
510	IF( lwp ) WRITE(numout, *) 'kt =', kt
511	IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
512	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
513	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
514	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
515	END IF
516	!
517	zht(:,:) = 0.0_wp
518	DO jk = 1, jpkm1
519	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
520	END DO
521	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
522	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
523	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
524	!
525	zht(:,:) = 0.0_wp
526	DO jk = 1, jpkm1
527	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
528	END DO
529	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
530	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
531	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
532	!
533	zht(:,:) = 0.0_wp
534	DO jk = 1, jpkm1
535	zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
536	END DO
537	z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
538	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
539	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
540	!
541	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) )
542	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
543	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
544	!
545	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) )
546	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
547	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
548	!
549	z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) )
550	CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain
551	IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
552	END IF
553
554	! *********************************** !
555	! After scale factors at u- v- points !
556	! *********************************** !
557
558	CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
559	CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
560
561	! *********************************** !
562	! After depths at u- v points !
563	! *********************************** !
564
565	hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
566	hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
567	DO jk = 2, jpkm1
568	hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
569	hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
570	END DO
571	! ! Inverse of the local depth
572	!!gm BUG ? don't understand the use of umask_i here .....
573	r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
574	r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
575	!
576	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt')
577	!
578	END SUBROUTINE dom_vvl_sf_nxt
579
580
581	SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
582	!!----------------------------------------------------------------------
583	!! * ROUTINE dom_vvl_sf_update *
584	!!
585	!! ** Purpose : for z tilde case: compute time filter and swap of scale factors
586	!! compute all depths and related variables for next time step
587	!! write outputs and restart file
588	!!
589	!! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
590	!! - reconstruct scale factor at other grid points (interpolate)
591	!! - recompute depths and water height fields
592	!!
593	!! ** Action : - tilde_e3t_(b/n) ready for next time step
594	!! - Recompute:
595	!! e3(u/v)_b
596	!! e3w(:,:,:,Kmm)
597	!! e3(u/v)w_b
598	!! e3(u/v)w_n
599	!! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
600	!! h(u/v) and h(u/v)r
601	!!
602	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
603	!! Leclair, M., and G. Madec, 2011, Ocean Modelling.
604	!!----------------------------------------------------------------------
605	INTEGER, INTENT( in ) :: kt ! time step
606	INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
607	!
608	INTEGER :: ji, jj, jk ! dummy loop indices
609	REAL(wp) :: zcoef ! local scalar
610	!!----------------------------------------------------------------------
611	!
612	IF( ln_linssh ) RETURN ! No calculation in linear free surface
613	!
614	IF( ln_timing ) CALL timing_start('dom_vvl_sf_update')
615	!
616	IF( kt == nit000 ) THEN
617	IF(lwp) WRITE(numout,*)
618	IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
619	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
620	ENDIF
621	!
622	! Time filter and swap of scale factors
623	! =====================================
624	! - ML - e3(t/u/v)_b are allready computed in dynnxt.
625	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
626	IF( neuler == 0 .AND. kt == nit000 ) THEN
627	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
628	ELSE
629	tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &
630	& + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
631	ENDIF
632	tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
633	ENDIF
634
635	! Compute all missing vertical scale factor and depths
636	! ====================================================
637	! Horizontal scale factor interpolations
638	! --------------------------------------
639	! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
640	! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
641
642	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
643
644	! Vertical scale factor interpolations
645	CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )
646	CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
647	CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
648	CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' )
649	CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
650	CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
651
652	! t- and w- points depth (set the isf depth as it is in the initial step)
653	gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
654	gdepw(:,:,1,Kmm) = 0.0_wp
655	gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
656	DO_3D_11_11( 2, jpk )
657	! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
658	! 1 for jk = mikt
659	zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
660	gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
661	gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) &
662	& + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) )
663	gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
664	END_3D
665
666	! Local depth and Inverse of the local depth of the water
667	! -------------------------------------------------------
668	!
669	ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
670	DO jk = 2, jpkm1
671	ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
672	END DO
673
674	! write restart file
675	! ==================
676	IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
677	!
678	IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update')
679	!
680	END SUBROUTINE dom_vvl_sf_update
681
682
683	SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
684	!!---------------------------------------------------------------------
685	!! * ROUTINE dom_vvl__interpol *
686	!!
687	!! ** Purpose : interpolate scale factors from one grid point to another
688	!!
689	!! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)
690	!! - horizontal interpolation: grid cell surface averaging
691	!! - vertical interpolation: simple averaging
692	!!----------------------------------------------------------------------
693	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated
694	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3
695	CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors
696	! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'
697	!
698	INTEGER :: ji, jj, jk ! dummy loop indices
699	REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F
700	!!----------------------------------------------------------------------
701	!
702	IF(ln_wd_il) THEN
703	zlnwd = 1.0_wp
704	ELSE
705	zlnwd = 0.0_wp
706	END IF
707	!
708	SELECT CASE ( pout ) !== type of interpolation ==!
709	!
710	CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean
711	DO_3D_10_10( 1, jpk )
712	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &
713	& * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &
714	& + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
715	END_3D
716	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
717	pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
718	!
719	CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean
720	DO_3D_10_10( 1, jpk )
721	pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &
722	& * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &
723	& + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
724	END_3D
725	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
726	pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
727	!
728	CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean
729	DO_3D_10_10( 1, jpk )
730	pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
731	& * r1_e1e2f(ji,jj) &
732	& * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &
733	& + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
734	END_3D
735	CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
736	pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
737	!
738	CASE( 'W' ) !* from T- to W-point : vertical simple mean
739	!
740	pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
741	! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
742	!!gm BUG? use here wmask in case of ISF ? to be checked
743	DO jk = 2, jpk
744	pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
745	& * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &
746	& + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
747	& * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )
748	END DO
749	!
750	CASE( 'UW' ) !* from U- to UW-point : vertical simple mean
751	!
752	pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
753	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
754	!!gm BUG? use here wumask in case of ISF ? to be checked
755	DO jk = 2, jpk
756	pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
757	& * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &
758	& + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
759	& * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )
760	END DO
761	!
762	CASE( 'VW' ) !* from V- to VW-point : vertical simple mean
763	!
764	pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
765	! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
766	!!gm BUG? use here wvmask in case of ISF ? to be checked
767	DO jk = 2, jpk
768	pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &
769	& * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &
770	& + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &
771	& * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )
772	END DO
773	END SELECT
774	!
775	END SUBROUTINE dom_vvl_interpol
776
777
778	SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
779	!!---------------------------------------------------------------------
780	!! * ROUTINE dom_vvl_rst *
781	!!
782	!! ** Purpose : Read or write VVL file in restart file
783	!!
784	!! ** Method : use of IOM library
785	!! if the restart does not contain vertical scale factors,
786	!! they are set to the _0 values
787	!! if the restart does not contain vertical scale factors increments (z_tilde),
788	!! they are set to 0.
789	!!----------------------------------------------------------------------
790	INTEGER , INTENT(in) :: kt ! ocean time-step
791	INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
792	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
793	!
794	INTEGER :: ji, jj, jk
795	INTEGER :: id1, id2, id3, id4, id5 ! local integers
796	!!----------------------------------------------------------------------
797	!
798	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
799	! ! ===============
800	IF( ln_rstart ) THEN !* Read the restart file
801	CALL rst_read_open ! open the restart file if necessary
802	CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios )
803	!
804	id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
805	id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
806	id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
807	id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
808	id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
809	!
810	! ! --------- !
811	! ! all cases !
812	! ! --------- !
813	!
814	IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist
815	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
816	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
817	! needed to restart if land processor not computed
818	IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
819	WHERE ( tmask(:,:,:) == 0.0_wp )
820	e3t(:,:,:,Kmm) = e3t_0(:,:,:)
821	e3t(:,:,:,Kbb) = e3t_0(:,:,:)
822	END WHERE
823	IF( neuler == 0 ) THEN
824	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
825	ENDIF
826	ELSE IF( id1 > 0 ) THEN
827	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
828	IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
829	IF(lwp) write(numout,*) 'neuler is forced to 0'
830	CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
831	e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
832	neuler = 0
833	ELSE IF( id2 > 0 ) THEN
834	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
835	IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
836	IF(lwp) write(numout,*) 'neuler is forced to 0'
837	CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
838	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
839	neuler = 0
840	ELSE
841	IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
842	IF(lwp) write(numout,*) 'Compute scale factor from sshn'
843	IF(lwp) write(numout,*) 'neuler is forced to 0'
844	DO jk = 1, jpk
845	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
846	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
847	& + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
848	END DO
849	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
850	neuler = 0
851	ENDIF
852	! ! ----------- !
853	IF( ln_vvl_zstar ) THEN ! z_star case !
854	! ! ----------- !
855	IF( MIN( id3, id4 ) > 0 ) THEN
856	CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
857	ENDIF
858	! ! ----------------------- !
859	ELSE ! z_tilde and layer cases !
860	! ! ----------------------- !
861	IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist
862	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
863	CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
864	ELSE ! one at least array is missing
865	tilde_e3t_b(:,:,:) = 0.0_wp
866	tilde_e3t_n(:,:,:) = 0.0_wp
867	ENDIF
868	! ! ------------ !
869	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
870	! ! ------------ !
871	IF( id5 > 0 ) THEN ! required array exists
872	CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
873	ELSE ! array is missing
874	hdiv_lf(:,:,:) = 0.0_wp
875	ENDIF
876	ENDIF
877	ENDIF
878	!
879	ELSE !* Initialize at "rest"
880	!
881
882	IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential
883	!
884	IF( cn_cfg == 'wad' ) THEN
885	! Wetting and drying test case
886	CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) )
887	ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones
888	ssh (:,:,Kmm) = ssh(:,:,Kbb)
889	uu (:,:,:,Kmm) = uu (:,:,:,Kbb)
890	vv (:,:,:,Kmm) = vv (:,:,:,Kbb)
891	ELSE
892	! if not test case
893	ssh(:,:,Kmm) = -ssh_ref
894	ssh(:,:,Kbb) = -ssh_ref
895
896	DO_2D_11_11
897	IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth
898	ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
899	ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
900	ENDIF
901	END_2D
902	ENDIF !If test case else
903
904	! Adjust vertical metrics for all wad
905	DO jk = 1, jpk
906	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
907	& / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
908	& + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
909	END DO
910	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
911
912	DO ji = 1, jpi
913	DO jj = 1, jpj
914	IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
915	CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
916	ENDIF
917	END DO
918	END DO
919	!
920	ELSE
921	!
922	! Just to read set ssh in fact, called latter once vertical grid
923	! is set up:
924	! CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) )
925	! !
926	! DO jk=1,jpk
927	! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
928	! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
929	! END DO
930	! e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
931	ssh(:,:,Kmm)=0._wp
932	e3t(:,:,:,Kmm)=e3t_0(:,:,:)
933	e3t(:,:,:,Kbb)=e3t_0(:,:,:)
934	!
935	END IF ! end of ll_wd edits
936
937	IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
938	tilde_e3t_b(:,:,:) = 0._wp
939	tilde_e3t_n(:,:,:) = 0._wp
940	IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
941	END IF
942	ENDIF
943	!
944	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
945	! ! ===================
946	IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
947	IF( lwxios ) CALL iom_swap( cwxios_context )
948	! ! --------- !
949	! ! all cases !
950	! ! --------- !
951	CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
952	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
953	! ! ----------------------- !
954	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
955	! ! ----------------------- !
956	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
957	CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
958	END IF
959	! ! -------------!
960	IF( ln_vvl_ztilde ) THEN ! z_tilde case !
961	! ! ------------ !
962	CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
963	ENDIF
964	!
965	IF( lwxios ) CALL iom_swap( cxios_context )
966	ENDIF
967	!
968	END SUBROUTINE dom_vvl_rst
969
970
971	SUBROUTINE dom_vvl_ctl
972	!!---------------------------------------------------------------------
973	!! * ROUTINE dom_vvl_ctl *
974	!!
975	!! ** Purpose : Control the consistency between namelist options
976	!! for vertical coordinate
977	!!----------------------------------------------------------------------
978	INTEGER :: ioptio, ios
979	!!
980	NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
981	& ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &
982	& rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe
983	!!----------------------------------------------------------------------
984	!
985	READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
986	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
987	READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
988	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
989	IF(lwm) WRITE ( numond, nam_vvl )
990	!
991	IF(lwp) THEN ! Namelist print
992	WRITE(numout,*)
993	WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
994	WRITE(numout,*) '~~~~~~~~~~~'
995	WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate'
996	WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar
997	WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde
998	WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer
999	WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar
1000	WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor
1001	WRITE(numout,*) ' !'
1002	WRITE(numout,*) ' thickness diffusion coefficient rn_ahe3 = ', rn_ahe3
1003	WRITE(numout,*) ' maximum e3t deformation fractional change rn_zdef_max = ', rn_zdef_max
1004	IF( ln_vvl_ztilde_as_zstar ) THEN
1005	WRITE(numout,*) ' ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1006	WRITE(numout,*) ' ignoring namelist timescale parameters and using:'
1007	WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'
1008	WRITE(numout,*) ' rn_rst_e3t = 0.e0'
1009	WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1010	WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt'
1011	ELSE
1012	WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t
1013	WRITE(numout,*) ' z-tilde cutoff frequency of low-pass filter (days) rn_lf_cutoff = ', rn_lf_cutoff
1014	ENDIF
1015	WRITE(numout,*) ' debug prints flag ln_vvl_dbg = ', ln_vvl_dbg
1016	ENDIF
1017	!
1018	ioptio = 0 ! Parameter control
1019	IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1020	IF( ln_vvl_zstar ) ioptio = ioptio + 1
1021	IF( ln_vvl_ztilde ) ioptio = ioptio + 1
1022	IF( ln_vvl_layer ) ioptio = ioptio + 1
1023	!
1024	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1025	!
1026	IF(lwp) THEN ! Print the choice
1027	WRITE(numout,*)
1028	IF( ln_vvl_zstar ) WRITE(numout,*) ' ==>>> zstar vertical coordinate is used'
1029	IF( ln_vvl_ztilde ) WRITE(numout,*) ' ==>>> ztilde vertical coordinate is used'
1030	IF( ln_vvl_layer ) WRITE(numout,*) ' ==>>> layer vertical coordinate is used'
1031	IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' ==>>> to emulate a zstar coordinate'
1032	ENDIF
1033	!
1034	#if defined key_agrif
1035	IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1036	#endif
1037	!
1038	END SUBROUTINE dom_vvl_ctl
1039
1040	!!======================================================================
1041	END MODULE domvvl

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/OCE/DOM/domvvl.F90 @ 12377

Download in other formats: