Context Navigation

domvvl.F90 @ 14585

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

Download in other formats: