Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

domvvl.F90

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

Download in other formats: