Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domvvl.F90 @ 7910

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

Download in other formats: