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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domzgr.F90 in branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5802

Last change on this file since 5802 was 5802, checked in by mathiot, 9 years ago
ice sheet coupling: reproducibility fixed in MISOMIP configuration + contrain bathy to be constant during all the run
Property svn:keywords set to `Id`
File size: 135.0 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
7	!! ! 1997-07 (G. Madec) lbc_lnk call
8	!! ! 1997-04 (J.-O. Beismann)
9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
19	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
20	!!----------------------------------------------------------------------
21
22	!!----------------------------------------------------------------------
23	!! dom_zgr : defined the ocean vertical coordinate system
24	!! zgr_bat : bathymetry fields (levels and meters)
25	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
26	!! zgr_bat_ctl : check the bathymetry files
27	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
28	!! zgr_z : reference z-coordinate
29	!! zgr_zco : z-coordinate
30	!! zgr_zps : z-coordinate with partial steps
31	!! zgr_sco : s-coordinate
32	!! fssig : tanh stretch function
33	!! fssig1 : Song and Haidvogel 1994 stretch function
34	!! fgamma : Siddorn and Furner 2012 stretching function
35	!!---------------------------------------------------------------------
36	USE oce ! ocean variables
37	USE dom_oce ! ocean domain
38	USE closea ! closed seas
39	USE c1d ! 1D vertical configuration
40	USE in_out_manager ! I/O manager
41	USE iom ! I/O library
42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
43	USE lib_mpp ! distributed memory computing library
44	USE wrk_nemo ! Memory allocation
45	USE timing ! Timing
46
47	IMPLICIT NONE
48	PRIVATE
49
50	PUBLIC dom_zgr ! called by dom_init.F90
51
52	! !!* Namelist namzgr_sco *
53	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
54	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
55	!
56	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
57	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
58	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
59	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
60	! Song and Haidvogel 1994 stretching parameters
61	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
62	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
63	REAL(wp) :: rn_bb ! stretching parameter
64	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
65	! Siddorn and Furner stretching parameters
66	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
67	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
68	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
69	REAL(wp) :: rn_zs ! depth of surface grid box
70	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
71	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
72	REAL(wp) :: rn_zb_b ! offset for calculating Zb
73
74	!! * Substitutions
75	# include "domzgr_substitute.h90"
76	# include "vectopt_loop_substitute.h90"
77	!!----------------------------------------------------------------------
78	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
79	!! $Id$
80	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
81	!!----------------------------------------------------------------------
82	CONTAINS
83
84	SUBROUTINE dom_zgr
85	!!----------------------------------------------------------------------
86	!! * ROUTINE dom_zgr *
87	!!
88	!! ** Purpose : set the depth of model levels and the resulting
89	!! vertical scale factors.
90	!!
91	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
92	!! - read/set ocean depth and ocean levels (bathy, mbathy)
93	!! - vertical coordinate (gdep., e3.) depending on the
94	!! coordinate chosen :
95	!! ln_zco=T z-coordinate
96	!! ln_zps=T z-coordinate with partial steps
97	!! ln_zco=T s-coordinate
98	!!
99	!! ** Action : define gdep., e3., mbathy and bathy
100	!!----------------------------------------------------------------------
101	INTEGER :: ioptio, ibat ! local integer
102	INTEGER :: ios
103	!
104	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
105	!!----------------------------------------------------------------------
106	!
107	IF( nn_timing == 1 ) CALL timing_start('dom_zgr')
108	!
109	REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
110	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
111	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
112
113	REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
114	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
115	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
116	IF(lwm) WRITE ( numond, namzgr )
117
118	IF(lwp) THEN ! Control print
119	WRITE(numout,*)
120	WRITE(numout,*) 'dom_zgr : vertical coordinate'
121	WRITE(numout,*) '~~~~~~~'
122	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
123	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
124	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
125	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
126	WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
127	ENDIF
128
129	ioptio = 0 ! Check Vertical coordinate options
130	IF( ln_zco ) ioptio = ioptio + 1
131	IF( ln_zps ) ioptio = ioptio + 1
132	IF( ln_sco ) ioptio = ioptio + 1
133	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
134	!
135	! Build the vertical coordinate system
136	! ------------------------------------
137	CALL zgr_z ! Reference z-coordinate system (always called)
138	CALL zgr_bat ! Bathymetry fields (levels and meters)
139	IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain
140	IF( ln_zco ) CALL zgr_zco ! z-coordinate
141	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
142	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
143	!
144	! final adjustment of mbathy & check
145	! -----------------------------------
146	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
147	IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
148	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
149	CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
150	!
151	IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain
152	ibat = mbathy(2,2)
153	mbathy(:,:) = ibat
154	END IF
155	!
156	IF( nprint == 1 .AND. lwp ) THEN
157	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
158	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
159	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
160	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
161	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
162	& ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)), &
163	& ' w ', MINVAL( e3w_0(:,:,:) )
164
165	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
166	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
167	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
168	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
169	& ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)), &
170	& ' w ', MAXVAL( e3w_0(:,:,:) )
171	ENDIF
172	!
173	IF( nn_timing == 1 ) CALL timing_stop('dom_zgr')
174	!
175	END SUBROUTINE dom_zgr
176
177
178	SUBROUTINE zgr_z
179	!!----------------------------------------------------------------------
180	!! * ROUTINE zgr_z *
181	!!
182	!! ** Purpose : set the depth of model levels and the resulting
183	!! vertical scale factors.
184	!!
185	!! ** Method : z-coordinate system (use in all type of coordinate)
186	!! The depth of model levels is defined from an analytical
187	!! function the derivative of which gives the scale factors.
188	!! both depth and scale factors only depend on k (1d arrays).
189	!! w-level: gdepw_1d = gdep(k)
190	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
191	!! t-level: gdept_1d = gdep(k+0.5)
192	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
193	!!
194	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
195	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
196	!!
197	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
198	!!----------------------------------------------------------------------
199	INTEGER :: jk ! dummy loop indices
200	REAL(wp) :: zt, zw ! temporary scalars
201	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
202	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
203	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
204	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
205	!!----------------------------------------------------------------------
206	!
207	IF( nn_timing == 1 ) CALL timing_start('zgr_z')
208	!
209	! Set variables from parameters
210	! ------------------------------
211	zkth = ppkth ; zacr = ppacr
212	zdzmin = ppdzmin ; zhmax = pphmax
213	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
214
215	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
216	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
217	IF( ppa1 == pp_to_be_computed .AND. &
218	& ppa0 == pp_to_be_computed .AND. &
219	& ppsur == pp_to_be_computed ) THEN
220	!
221	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
222	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
223	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
224	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
225	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
226	ELSE
227	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
228	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
229	ENDIF
230
231	IF(lwp) THEN ! Parameter print
232	WRITE(numout,*)
233	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
234	WRITE(numout,*) ' ~~~~~~~'
235	IF( ppkth == 0._wp ) THEN
236	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
237	WRITE(numout,*) ' Total depth :', zhmax
238	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
239	ELSE
240	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
241	WRITE(numout,*) ' zsur, za0, za1 computed from '
242	WRITE(numout,*) ' zdzmin = ', zdzmin
243	WRITE(numout,*) ' zhmax = ', zhmax
244	ENDIF
245	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
246	WRITE(numout,*) ' zsur = ', zsur
247	WRITE(numout,*) ' za0 = ', za0
248	WRITE(numout,*) ' za1 = ', za1
249	WRITE(numout,*) ' zkth = ', zkth
250	WRITE(numout,*) ' zacr = ', zacr
251	IF( ldbletanh ) THEN
252	WRITE(numout,*) ' (Double tanh za2 = ', za2
253	WRITE(numout,*) ' parameters) zkth2= ', zkth2
254	WRITE(numout,*) ' zacr2= ', zacr2
255	ENDIF
256	ENDIF
257	ENDIF
258
259
260	! Reference z-coordinate (depth - scale factor at T- and W-points)
261	! ======================
262	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
263	za1 = zhmax / FLOAT(jpk-1)
264	DO jk = 1, jpk
265	zw = FLOAT( jk )
266	zt = FLOAT( jk ) + 0.5_wp
267	gdepw_1d(jk) = ( zw - 1 ) * za1
268	gdept_1d(jk) = ( zt - 1 ) * za1
269	e3w_1d (jk) = za1
270	e3t_1d (jk) = za1
271	END DO
272	ELSE ! Madec & Imbard 1996 function
273	IF( .NOT. ldbletanh ) THEN
274	DO jk = 1, jpk
275	zw = REAL( jk , wp )
276	zt = REAL( jk , wp ) + 0.5_wp
277	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
278	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
279	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
280	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
281	END DO
282	ELSE
283	DO jk = 1, jpk
284	zw = FLOAT( jk )
285	zt = FLOAT( jk ) + 0.5_wp
286	! Double tanh function
287	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
288	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
289	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
290	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
291	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
292	& + za2 * TANH( (zw-zkth2) / zacr2 )
293	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
294	& + za2 * TANH( (zt-zkth2) / zacr2 )
295	END DO
296	ENDIF
297	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
298	ENDIF
299
300	IF ( ln_isfcav ) THEN
301	! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
302	! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
303	DO jk = 1, jpkm1
304	e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
305	END DO
306	e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
307
308	DO jk = 2, jpk
309	e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
310	END DO
311	e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
312	END IF
313
314	!!gm BUG in s-coordinate this does not work!
315	! deepest/shallowest W level Above/Below ~10m
316	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
317	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
318	nla10 = nlb10 - 1 ! deepest W level Above ~10m
319	!!gm end bug
320
321	IF(lwp) THEN ! control print
322	WRITE(numout,*)
323	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
324	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
325	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
326	ENDIF
327	DO jk = 1, jpk ! control positivity
328	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
329	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
330	END DO
331	!
332	IF( nn_timing == 1 ) CALL timing_stop('zgr_z')
333	!
334	END SUBROUTINE zgr_z
335
336
337	SUBROUTINE zgr_bat
338	!!----------------------------------------------------------------------
339	!! * ROUTINE zgr_bat *
340	!!
341	!! ** Purpose : set bathymetry both in levels and meters
342	!!
343	!! ** Method : read or define mbathy and bathy arrays
344	!! * level bathymetry:
345	!! The ocean basin geometry is given by a two-dimensional array,
346	!! mbathy, which is defined as follow :
347	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
348	!! at t-point (ji,jj).
349	!! = 0 over the continental t-point.
350	!! The array mbathy is checked to verified its consistency with
351	!! model option. in particular:
352	!! mbathy must have at least 1 land grid-points (mbathy<=0)
353	!! along closed boundary.
354	!! mbathy must be cyclic IF jperio=1.
355	!! mbathy must be lower or equal to jpk-1.
356	!! isolated ocean grid points are suppressed from mbathy
357	!! since they are only connected to remaining
358	!! ocean through vertical diffusion.
359	!! ntopo=-1 : rectangular channel or bassin with a bump
360	!! ntopo= 0 : flat rectangular channel or basin
361	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
362	!! bathy is read in 'bathy_meter.nc' NetCDF file
363	!!
364	!! ** Action : - mbathy: level bathymetry (in level index)
365	!! - bathy : meter bathymetry (in meters)
366	!!----------------------------------------------------------------------
367	INTEGER :: ji, jj, jk ! dummy loop indices
368	INTEGER :: inum ! temporary logical unit
369	INTEGER :: ierror ! error flag
370	INTEGER :: ii_bump, ij_bump, ih ! bump center position
371	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
372	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
373	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
374	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data
375	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data
376	!!----------------------------------------------------------------------
377	!
378	IF( nn_timing == 1 ) CALL timing_start('zgr_bat')
379	!
380	IF(lwp) WRITE(numout,*)
381	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
382	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
383	! ! ================== !
384	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
385	! ! ================== !
386	! ! global domain level and meter bathymetry (idta,zdta)
387	!
388	ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
389	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
390	ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
391	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
392	!
393	IF( ntopo == 0 ) THEN ! flat basin
394	IF(lwp) WRITE(numout,*)
395	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
396	IF( rn_bathy > 0.01 ) THEN
397	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
398	zdta(:,:) = rn_bathy
399	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
400	idta(:,:) = jpkm1
401	ELSE ! z-coordinate (zco or zps): step-like topography
402	idta(:,:) = jpkm1
403	DO jk = 1, jpkm1
404	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
405	END DO
406	ENDIF
407	ELSE
408	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
409	idta(:,:) = jpkm1 ! before last level
410	zdta(:,:) = gdepw_1d(jpk) ! last w-point depth
411	h_oce = gdepw_1d(jpk)
412	ENDIF
413	ELSE ! bump centered in the basin
414	IF(lwp) WRITE(numout,*)
415	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
416	ii_bump = jpidta / 2 ! i-index of the bump center
417	ij_bump = jpjdta / 2 ! j-index of the bump center
418	r_bump = 50000._wp ! bump radius (meters)
419	h_bump = 2700._wp ! bump height (meters)
420	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
421	IF(lwp) WRITE(numout,*) ' bump characteristics: '
422	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
423	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
424	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
425	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
426	!
427	DO jj = 1, jpjdta ! zdta :
428	DO ji = 1, jpidta
429	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
430	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
431	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
432	END DO
433	END DO
434	! ! idta :
435	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
436	idta(:,:) = jpkm1
437	ELSE ! z-coordinate (zco or zps): step-like topography
438	idta(:,:) = jpkm1
439	DO jk = 1, jpkm1
440	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
441	END DO
442	ENDIF
443	ENDIF
444	! ! set GLOBAL boundary conditions
445	! ! Caution : idta on the global domain: use of jperio, not nperio
446	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
447	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
448	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
449	ELSEIF( jperio == 2 ) THEN
450	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
451	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
452	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
453	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp
454	ELSE
455	ih = 0 ; zh = 0._wp
456	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
457	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
458	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
459	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
460	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
461	ENDIF
462
463	! ! local domain level and meter bathymetries (mbathy,bathy)
464	mbathy(:,:) = 0 ! set to zero extra halo points
465	bathy (:,:) = 0._wp ! (require for mpp case)
466	DO jj = 1, nlcj ! interior values
467	DO ji = 1, nlci
468	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
469	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
470	END DO
471	END DO
472	risfdep(:,:)=0.e0
473	misfdep(:,:)=1
474	!
475	DEALLOCATE( idta, zdta )
476	!
477	! ! ================ !
478	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
479	! ! ================ !
480	!
481	IF( ln_zco ) THEN ! zco : read level bathymetry
482	CALL iom_open ( 'bathy_level.nc', inum )
483	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
484	CALL iom_close( inum )
485	mbathy(:,:) = INT( bathy(:,:) )
486	! ! =====================
487	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
488	! ! =====================
489	IF( nn_cla == 0 ) THEN
490	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
491	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
492	DO ji = mi0(ii0), mi1(ii1)
493	DO jj = mj0(ij0), mj1(ij1)
494	mbathy(ji,jj) = 15
495	END DO
496	END DO
497	IF(lwp) WRITE(numout,*)
498	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
499	!
500	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
501	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
502	DO ji = mi0(ii0), mi1(ii1)
503	DO jj = mj0(ij0), mj1(ij1)
504	mbathy(ji,jj) = 12
505	END DO
506	END DO
507	IF(lwp) WRITE(numout,*)
508	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
509	ENDIF
510	!
511	ENDIF
512	!
513	ENDIF
514	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
515	CALL iom_open ( 'bathy_meter.nc', inum )
516	IF ( ln_isfcav ) THEN
517	CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
518	ELSE
519	CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr )
520	END IF
521	CALL iom_close( inum )
522	!
523	risfdep(:,:)=0._wp
524	misfdep(:,:)=1
525	IF ( ln_isfcav ) THEN
526	CALL iom_open ( 'isf_draft_meter.nc', inum )
527	CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep )
528	CALL iom_close( inum )
529	WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
530
531	! set grounded point to 0 (treshold at 1cm, have to be update after first coupling experience)
532	WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 )
533	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
534	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
535	END WHERE
536	END IF
537	!
538	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
539	!
540	IF( nn_cla == 0 ) THEN
541	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
542	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
543	DO ji = mi0(ii0), mi1(ii1)
544	DO jj = mj0(ij0), mj1(ij1)
545	bathy(ji,jj) = 284._wp
546	END DO
547	END DO
548	IF(lwp) WRITE(numout,*)
549	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
550	!
551	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
552	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
553	DO ji = mi0(ii0), mi1(ii1)
554	DO jj = mj0(ij0), mj1(ij1)
555	bathy(ji,jj) = 137._wp
556	END DO
557	END DO
558	IF(lwp) WRITE(numout,*)
559	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
560	ENDIF
561	!
562	ENDIF
563	!
564	ENDIF
565	! ! =============== !
566	ELSE ! error !
567	! ! =============== !
568	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
569	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
570	ENDIF
571	!
572	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
573	!
574	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
575	! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
576	IF ( ln_isfcav ) THEN
577	WHERE (bathy == risfdep)
578	bathy = 0.0_wp ; risfdep = 0.0_wp
579	END WHERE
580	END IF
581	! end patch
582	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
583	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
584	ENDIF
585	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
586	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
587	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
588	END WHERE
589	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
590	ENDIF
591	!
592	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')
593	!
594	END SUBROUTINE zgr_bat
595
596
597	SUBROUTINE zgr_bat_zoom
598	!!----------------------------------------------------------------------
599	!! * ROUTINE zgr_bat_zoom *
600	!!
601	!! ** Purpose : - Close zoom domain boundary if necessary
602	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
603	!!
604	!! ** Method :
605	!!
606	!! ** Action : - update mbathy: level bathymetry (in level index)
607	!!----------------------------------------------------------------------
608	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
609	!!----------------------------------------------------------------------
610	!
611	IF(lwp) WRITE(numout,*)
612	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
613	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
614	!
615	! Zoom domain
616	! ===========
617	!
618	! Forced closed boundary if required
619	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
620	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
621	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
622	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
623	!
624	! Configuration specific domain modifications
625	! (here, ORCA arctic configuration: suppress Med Sea)
626	IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
627	SELECT CASE ( jp_cfg )
628	! ! =======================
629	CASE ( 2 ) ! ORCA_R2 configuration
630	! ! =======================
631	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
632	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
633	ij0 = 98 ; ij1 = 110
634	! ! =======================
635	CASE ( 05 ) ! ORCA_R05 configuration
636	! ! =======================
637	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
638	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
639	ij0 = 314 ; ij1 = 370
640	END SELECT
641	!
642	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
643	!
644	ENDIF
645	!
646	END SUBROUTINE zgr_bat_zoom
647
648
649	SUBROUTINE zgr_bat_ctl
650	!!----------------------------------------------------------------------
651	!! * ROUTINE zgr_bat_ctl *
652	!!
653	!! ** Purpose : check the bathymetry in levels
654	!!
655	!! ** Method : The array mbathy is checked to verified its consistency
656	!! with the model options. in particular:
657	!! mbathy must have at least 1 land grid-points (mbathy<=0)
658	!! along closed boundary.
659	!! mbathy must be cyclic IF jperio=1.
660	!! mbathy must be lower or equal to jpk-1.
661	!! isolated ocean grid points are suppressed from mbathy
662	!! since they are only connected to remaining
663	!! ocean through vertical diffusion.
664	!! C A U T I O N : mbathy will be modified during the initializa-
665	!! tion phase to become the number of non-zero w-levels of a water
666	!! column, with a minimum value of 1.
667	!!
668	!! ** Action : - update mbathy: level bathymetry (in level index)
669	!! - update bathy : meter bathymetry (in meters)
670	!!----------------------------------------------------------------------
671	!!
672	INTEGER :: ji, jj, jl ! dummy loop indices
673	INTEGER :: icompt, ibtest, ikmax ! temporary integers
674	REAL(wp), POINTER, DIMENSION(:,:) :: zbathy
675
676	!!----------------------------------------------------------------------
677	!
678	IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl')
679	!
680	CALL wrk_alloc( jpi, jpj, zbathy )
681	!
682	IF(lwp) WRITE(numout,*)
683	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
684	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
685	! ! Suppress isolated ocean grid points
686	IF(lwp) WRITE(numout,*)
687	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
688	IF(lwp) WRITE(numout,*)' -----------------------------------'
689	icompt = 0
690	DO jl = 1, 2
691	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
692	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
693	mbathy(jpi,:) = mbathy( 2 ,:)
694	ENDIF
695	DO jj = 2, jpjm1
696	DO ji = 2, jpim1
697	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
698	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
699	IF( ibtest < mbathy(ji,jj) ) THEN
700	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
701	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
702	mbathy(ji,jj) = ibtest
703	icompt = icompt + 1
704	ENDIF
705	END DO
706	END DO
707	END DO
708	IF( lk_mpp ) CALL mpp_sum( icompt )
709	IF( icompt == 0 ) THEN
710	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
711	ELSE
712	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
713	ENDIF
714	IF( lk_mpp ) THEN
715	zbathy(:,:) = FLOAT( mbathy(:,:) )
716	CALL lbc_lnk( zbathy, 'T', 1._wp )
717	mbathy(:,:) = INT( zbathy(:,:) )
718	ENDIF
719	! ! East-west cyclic boundary conditions
720	IF( nperio == 0 ) THEN
721	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
722	IF( lk_mpp ) THEN
723	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
724	IF( jperio /= 1 ) mbathy(1,:) = 0
725	ENDIF
726	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
727	IF( jperio /= 1 ) mbathy(nlci,:) = 0
728	ENDIF
729	ELSE
730	IF( ln_zco .OR. ln_zps ) THEN
731	mbathy( 1 ,:) = 0
732	mbathy(jpi,:) = 0
733	ELSE
734	mbathy( 1 ,:) = jpkm1
735	mbathy(jpi,:) = jpkm1
736	ENDIF
737	ENDIF
738	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
739	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
740	mbathy( 1 ,:) = mbathy(jpim1,:)
741	mbathy(jpi,:) = mbathy( 2 ,:)
742	ELSEIF( nperio == 2 ) THEN
743	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
744	ELSE
745	IF(lwp) WRITE(numout,*) ' e r r o r'
746	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
747	! STOP 'dom_mba'
748	ENDIF
749	! Boundary condition on mbathy
750	IF( .NOT.lk_mpp ) THEN
751	!!gm !!bug ??? think about it !
752	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
753	zbathy(:,:) = FLOAT( mbathy(:,:) )
754	CALL lbc_lnk( zbathy, 'T', 1._wp )
755	mbathy(:,:) = INT( zbathy(:,:) )
756	ENDIF
757	! Number of ocean level inferior or equal to jpkm1
758	ikmax = 0
759	DO jj = 1, jpj
760	DO ji = 1, jpi
761	ikmax = MAX( ikmax, mbathy(ji,jj) )
762	END DO
763	END DO
764	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
765	IF( ikmax > jpkm1 ) THEN
766	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
767	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
768	ELSE IF( ikmax < jpkm1 ) THEN
769	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
770	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
771	ENDIF
772
773	IF( lwp .AND. nprint == 1 ) THEN ! control print
774	WRITE(numout,*)
775	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
776	WRITE(numout,*) ' ------------------'
777	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
778	WRITE(numout,*)
779	ENDIF
780	!
781	CALL wrk_dealloc( jpi, jpj, zbathy )
782	!
783	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl')
784	!
785	END SUBROUTINE zgr_bat_ctl
786
787
788	SUBROUTINE zgr_bot_level
789	!!----------------------------------------------------------------------
790	!! * ROUTINE zgr_bot_level *
791	!!
792	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
793	!!
794	!! ** Method : computes from mbathy with a minimum value of 1 over land
795	!!
796	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
797	!! ocean level at t-, u- & v-points
798	!! (min value = 1 over land)
799	!!----------------------------------------------------------------------
800	!!
801	INTEGER :: ji, jj ! dummy loop indices
802	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
803	!!----------------------------------------------------------------------
804	!
805	IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level')
806	!
807	CALL wrk_alloc( jpi, jpj, zmbk )
808	!
809	IF(lwp) WRITE(numout,*)
810	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
811	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
812	!
813	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
814
815	! ! bottom k-index of W-level = mbkt+1
816	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
817	DO ji = 1, jpim1
818	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
819	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
820	END DO
821	END DO
822	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
823	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
824	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
825	!
826	CALL wrk_dealloc( jpi, jpj, zmbk )
827	!
828	IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level')
829	!
830	END SUBROUTINE zgr_bot_level
831
832	SUBROUTINE zgr_top_level
833	!!----------------------------------------------------------------------
834	!! * ROUTINE zgr_bot_level *
835	!!
836	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
837	!!
838	!! ** Method : computes from misfdep with a minimum value of 1
839	!!
840	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
841	!! ocean level at t-, u- & v-points
842	!! (min value = 1)
843	!!----------------------------------------------------------------------
844	!!
845	INTEGER :: ji, jj ! dummy loop indices
846	REAL(wp), POINTER, DIMENSION(:,:) :: zmik
847	!!----------------------------------------------------------------------
848	!
849	IF( nn_timing == 1 ) CALL timing_start('zgr_top_level')
850	!
851	CALL wrk_alloc( jpi, jpj, zmik )
852	!
853	IF(lwp) WRITE(numout,*)
854	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
855	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
856	!
857	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
858	! ! top k-index of W-level (=mikt)
859	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
860	DO ji = 1, jpim1
861	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
862	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
863	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
864	END DO
865	END DO
866
867	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
868	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
869	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
870	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
871	!
872	CALL wrk_dealloc( jpi, jpj, zmik )
873	!
874	IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level')
875	!
876	END SUBROUTINE zgr_top_level
877
878	SUBROUTINE zgr_zco
879	!!----------------------------------------------------------------------
880	!! * ROUTINE zgr_zco *
881	!!
882	!! ** Purpose : define the z-coordinate system
883	!!
884	!! ** Method : set 3D coord. arrays to reference 1D array
885	!!----------------------------------------------------------------------
886	INTEGER :: jk
887	!!----------------------------------------------------------------------
888	!
889	IF( nn_timing == 1 ) CALL timing_start('zgr_zco')
890	!
891	DO jk = 1, jpk
892	gdept_0 (:,:,jk) = gdept_1d(jk)
893	gdepw_0 (:,:,jk) = gdepw_1d(jk)
894	gdep3w_0(:,:,jk) = gdepw_1d(jk)
895	e3t_0 (:,:,jk) = e3t_1d (jk)
896	e3u_0 (:,:,jk) = e3t_1d (jk)
897	e3v_0 (:,:,jk) = e3t_1d (jk)
898	e3f_0 (:,:,jk) = e3t_1d (jk)
899	e3w_0 (:,:,jk) = e3w_1d (jk)
900	e3uw_0 (:,:,jk) = e3w_1d (jk)
901	e3vw_0 (:,:,jk) = e3w_1d (jk)
902	END DO
903	!
904	IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')
905	!
906	END SUBROUTINE zgr_zco
907
908
909	SUBROUTINE zgr_zps
910	!!----------------------------------------------------------------------
911	!! * ROUTINE zgr_zps *
912	!!
913	!! ** Purpose : the depth and vertical scale factor in partial step
914	!! z-coordinate case
915	!!
916	!! ** Method : Partial steps : computes the 3D vertical scale factors
917	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
918	!! a partial step representation of bottom topography.
919	!!
920	!! The reference depth of model levels is defined from an analytical
921	!! function the derivative of which gives the reference vertical
922	!! scale factors.
923	!! From depth and scale factors reference, we compute there new value
924	!! with partial steps on 3d arrays ( i, j, k ).
925	!!
926	!! w-level: gdepw_0(i,j,k) = gdep(k)
927	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
928	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
929	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
930	!!
931	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
932	!! we find the mbathy index of the depth at each grid point.
933	!! This leads us to three cases:
934	!!
935	!! - bathy = 0 => mbathy = 0
936	!! - 1 < mbathy < jpkm1
937	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
938	!!
939	!! Then, for each case, we find the new depth at t- and w- levels
940	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
941	!! and f-points.
942	!!
943	!! This routine is given as an example, it must be modified
944	!! following the user s desiderata. nevertheless, the output as
945	!! well as the way to compute the model levels and scale factors
946	!! must be respected in order to insure second order accuracy
947	!! schemes.
948	!!
949	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
950	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
951	!!
952	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
953	!!----------------------------------------------------------------------
954	!!
955	INTEGER :: ji, jj, jk ! dummy loop indices
956	INTEGER :: ik, it, ikb, ikt ! temporary integers
957	LOGICAL :: ll_print ! Allow control print for debugging
958	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
959	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
960	REAL(wp) :: zdiff ! temporary scalar
961	REAL(wp) :: zmax ! temporary scalar
962	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
963	!!---------------------------------------------------------------------
964	!
965	IF( nn_timing == 1 ) CALL timing_start('zgr_zps')
966	!
967	CALL wrk_alloc( jpi, jpj, jpk, zprt )
968	!
969	IF(lwp) WRITE(numout,*)
970	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
971	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
972	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
973
974	ll_print = .FALSE. ! Local variable for debugging
975
976	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
977	WRITE(numout,*)
978	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
979	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
980	ENDIF
981
982
983	! bathymetry in level (from bathy_meter)
984	! ===================
985	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
986	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
987	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
988	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
989	END WHERE
990
991	! Compute mbathy for ocean points (i.e. the number of ocean levels)
992	! find the number of ocean levels such that the last level thickness
993	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
994	! e3t_1d is the reference level thickness
995	DO jk = jpkm1, 1, -1
996	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
997	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
998	END DO
999
1000	! Scale factors and depth at T- and W-points
1001	DO jk = 1, jpk ! intitialization to the reference z-coordinate
1002	gdept_0(:,:,jk) = gdept_1d(jk)
1003	gdepw_0(:,:,jk) = gdepw_1d(jk)
1004	e3t_0 (:,:,jk) = e3t_1d (jk)
1005	e3w_0 (:,:,jk) = e3w_1d (jk)
1006	END DO
1007
1008	! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
1009	IF ( ln_isfcav ) CALL zgr_isf
1010
1011	! Scale factors and depth at T- and W-points
1012	IF ( .NOT. ln_isfcav ) THEN
1013	DO jj = 1, jpj
1014	DO ji = 1, jpi
1015	ik = mbathy(ji,jj)
1016	IF( ik > 0 ) THEN ! ocean point only
1017	! max ocean level case
1018	IF( ik == jpkm1 ) THEN
1019	zdepwp = bathy(ji,jj)
1020	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1021	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1022	e3t_0(ji,jj,ik ) = ze3tp
1023	e3t_0(ji,jj,ik+1) = ze3tp
1024	e3w_0(ji,jj,ik ) = ze3wp
1025	e3w_0(ji,jj,ik+1) = ze3tp
1026	gdepw_0(ji,jj,ik+1) = zdepwp
1027	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1028	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1029	!
1030	ELSE ! standard case
1031	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1032	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1033	ENDIF
1034	!gm Bug? check the gdepw_1d
1035	! ... on ik
1036	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1037	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1038	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1039	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1040	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1041	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1042	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1043	! ... on ik+1
1044	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1045	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1046	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1047	ENDIF
1048	ENDIF
1049	END DO
1050	END DO
1051	!
1052	it = 0
1053	DO jj = 1, jpj
1054	DO ji = 1, jpi
1055	ik = mbathy(ji,jj)
1056	IF( ik > 0 ) THEN ! ocean point only
1057	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1058	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1059	! test
1060	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1061	IF( zdiff <= 0._wp .AND. lwp ) THEN
1062	it = it + 1
1063	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1064	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1065	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1066	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1067	ENDIF
1068	ENDIF
1069	END DO
1070	END DO
1071	END IF
1072	!
1073	! Scale factors and depth at U-, V-, UW and VW-points
1074	DO jk = 1, jpk ! initialisation to z-scale factors
1075	e3u_0 (:,:,jk) = e3t_1d(jk)
1076	e3v_0 (:,:,jk) = e3t_1d(jk)
1077	e3uw_0(:,:,jk) = e3w_1d(jk)
1078	e3vw_0(:,:,jk) = e3w_1d(jk)
1079	END DO
1080
1081	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1082	DO jj = 1, jpjm1
1083	DO ji = 1, fs_jpim1 ! vector opt.
1084	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1085	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1086	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1087	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1088	END DO
1089	END DO
1090	END DO
1091	IF ( ln_isfcav ) THEN
1092	! (ISF) define e3uw (adapted for 2 cells in the water column)
1093	DO jj = 2, jpjm1
1094	DO ji = 2, fs_jpim1 ! vector opt.
1095	ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1096	ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1097	IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &
1098	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )
1099	ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1100	ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1101	IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &
1102	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )
1103	END DO
1104	END DO
1105	END IF
1106
1107	CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1108	CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1109	!
1110
1111	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1112	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1113	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1114	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1115	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1116	END DO
1117
1118	! Scale factor at F-point
1119	DO jk = 1, jpk ! initialisation to z-scale factors
1120	e3f_0(:,:,jk) = e3t_1d(jk)
1121	END DO
1122	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1123	DO jj = 1, jpjm1
1124	DO ji = 1, fs_jpim1 ! vector opt.
1125	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1126	END DO
1127	END DO
1128	END DO
1129	CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1130	!
1131	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1132	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1133	END DO
1134	!!gm bug ? : must be a do loop with mj0,mj1
1135	!
1136	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1137	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1138	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1139	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1140	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1141
1142	! Control of the sign
1143	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1144	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1145	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1146	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1147
1148	! Compute gdep3w_0 (vertical sum of e3w)
1149	IF ( ln_isfcav ) THEN ! if cavity
1150	WHERE (misfdep == 0) misfdep = 1
1151	DO jj = 1,jpj
1152	DO ji = 1,jpi
1153	gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1154	DO jk = 2, misfdep(ji,jj)
1155	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1156	END DO
1157	IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
1158	DO jk = misfdep(ji,jj) + 1, jpk
1159	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1160	END DO
1161	END DO
1162	END DO
1163	ELSE ! no cavity
1164	gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1165	DO jk = 2, jpk
1166	gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1167	END DO
1168	END IF
1169	! ! ================= !
1170	IF(lwp .AND. ll_print) THEN ! Control print !
1171	! ! ================= !
1172	DO jj = 1,jpj
1173	DO ji = 1, jpi
1174	ik = MAX( mbathy(ji,jj), 1 )
1175	zprt(ji,jj,1) = e3t_0 (ji,jj,ik)
1176	zprt(ji,jj,2) = e3w_0 (ji,jj,ik)
1177	zprt(ji,jj,3) = e3u_0 (ji,jj,ik)
1178	zprt(ji,jj,4) = e3v_0 (ji,jj,ik)
1179	zprt(ji,jj,5) = e3f_0 (ji,jj,ik)
1180	zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1181	END DO
1182	END DO
1183	WRITE(numout,*)
1184	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1185	WRITE(numout,*)
1186	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1187	WRITE(numout,*)
1188	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1189	WRITE(numout,*)
1190	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1191	WRITE(numout,*)
1192	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1193	WRITE(numout,*)
1194	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1195	ENDIF
1196	!
1197	CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1198	!
1199	IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
1200	!
1201	END SUBROUTINE zgr_zps
1202
1203	SUBROUTINE zgr_isf
1204	!!----------------------------------------------------------------------
1205	!! * ROUTINE zgr_isf *
1206	!!
1207	!! ** Purpose : check the bathymetry in levels
1208	!!
1209	!! ** Method : THe water column have to contained at least 2 cells
1210	!! Bathymetry and isfdraft are modified (dig/close) to respect
1211	!! this criterion.
1212	!!
1213	!!
1214	!! ** Action : - test compatibility between isfdraft and bathy
1215	!! - bathy and isfdraft are modified
1216	!!----------------------------------------------------------------------
1217	!!
1218	INTEGER :: ji, jj, jl, jk ! dummy loop indices
1219	INTEGER :: ik, it ! temporary integers
1220	INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
1221	REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t
1222	REAL(wp) :: zmax ! Maximum and minimum depth
1223	REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar
1224	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1225	REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t
1226	REAL(wp) :: zdiff ! temporary scalar
1227	REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)
1228	INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)
1229	!!---------------------------------------------------------------------
1230	!
1231	IF( nn_timing == 1 ) CALL timing_start('zgr_isf')
1232	!
1233	CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1234	CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1235
1236
1237	! (ISF) compute misfdep
1238	WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1
1239	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
1240	END WHERE
1241
1242	! Compute misfdep for ocean points (i.e. first wet level)
1243	! find the first ocean level such that the first level thickness
1244	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1245	! e3t_0 is the reference level thickness
1246	DO jk = 2, jpkm1
1247	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1248	WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
1249	END DO
1250	WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
1251	risfdep(:,:) = 0. ; misfdep(:,:) = 1
1252	END WHERE
1253
1254	! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
1255	icompt = 0
1256	! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1257	DO jl = 1, 10
1258	! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
1259	WHERE (bathy(:,:) .LE. risfdep(:,:)+1e-2 )
1260	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1261	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
1262	END WHERE
1263	WHERE (mbathy(:,:) <= 0)
1264	misfdep(:,:) = 0; risfdep(:,:) = 0._wp
1265	mbathy (:,:) = 0; bathy (:,:) = 0._wp
1266	ENDWHERE
1267	IF( lk_mpp ) THEN
1268	zbathy(:,:) = FLOAT( misfdep(:,:) )
1269	CALL lbc_lnk( zbathy, 'T', 1. )
1270	misfdep(:,:) = INT( zbathy(:,:) )
1271	CALL lbc_lnk( risfdep, 'T', 1. )
1272	CALL lbc_lnk( bathy, 'T', 1. )
1273	zbathy(:,:) = FLOAT( mbathy(:,:) )
1274	CALL lbc_lnk( zbathy, 'T', 1. )
1275	mbathy(:,:) = INT( zbathy(:,:) )
1276	ENDIF
1277	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1278	misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west
1279	misfdep(jpi,:) = misfdep( 2 ,:)
1280	ENDIF
1281
1282	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1283	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
1284	mbathy(jpi,:) = mbathy( 2 ,:)
1285	ENDIF
1286
1287	! split last cell if possible (only where water column is 2 cell or less)
1288	! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
1289	IF ( .NOT. ln_iscpl) THEN
1290	DO jk = jpkm1, 1, -1
1291	zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1292	WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1293	mbathy(:,:) = jk
1294	bathy(:,:) = zmax
1295	END WHERE
1296	END DO
1297	END IF
1298
1299	! split top cell if possible (only where water column is 2 cell or less)
1300	DO jk = 2, jpkm1
1301	zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1302	WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1303	misfdep(:,:) = jk
1304	risfdep(:,:) = zmax
1305	END WHERE
1306	END DO
1307
1308
1309	! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1310	DO jj = 1, jpj
1311	DO ji = 1, jpi
1312	! find the minimum change option:
1313	! test bathy
1314	IF (risfdep(ji,jj) .GT. 1) THEN
1315	IF ( .NOT. ln_iscpl ) THEN
1316	zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
1317	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1318	zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &
1319	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1320
1321	IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1322	IF (zbathydiff .LE. zrisfdepdiff) THEN
1323	bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1324	mbathy(ji,jj)= mbathy(ji,jj) + 1
1325	ELSE
1326	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1327	misfdep(ji,jj) = misfdep(ji,jj) - 1
1328	END IF
1329	ENDIF
1330	ELSE
1331	IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1332	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1333	misfdep(ji,jj) = misfdep(ji,jj) - 1
1334	END IF
1335	END IF
1336	END IF
1337	END DO
1338	END DO
1339
1340	! At least 2 levels for water thickness at T, U, and V point.
1341	DO jj = 1, jpj
1342	DO ji = 1, jpi
1343	! find the minimum change option:
1344	! test bathy
1345	IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1346	IF ( .NOT. ln_iscpl ) THEN
1347	zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
1348	& + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1349	zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) &
1350	& - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1351	IF (zbathydiff .LE. zrisfdepdiff) THEN
1352	mbathy(ji,jj) = mbathy(ji,jj) + 1
1353	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1354	ELSE
1355	misfdep(ji,jj)= misfdep(ji,jj) - 1
1356	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1357	END IF
1358	ELSE
1359	misfdep(ji,jj)= misfdep(ji,jj) - 1
1360	risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1361	END IF
1362	ENDIF
1363	END DO
1364	END DO
1365
1366	! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1367	DO jj = 1, jpjm1
1368	DO ji = 1, jpim1
1369	IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1370	IF ( .NOT. ln_iscpl ) THEN
1371	zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) &
1372	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))
1373	zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
1374	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1375	IF (zbathydiff .LE. zrisfdepdiff) THEN
1376	mbathy(ji,jj) = mbathy(ji,jj) + 1
1377	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )
1378	ELSE
1379	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
1380	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1381	END IF
1382	ELSE
1383	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
1384	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1385	END IF
1386	ENDIF
1387	END DO
1388	END DO
1389
1390	IF( lk_mpp ) THEN
1391	zbathy(:,:) = FLOAT( misfdep(:,:) )
1392	CALL lbc_lnk( zbathy, 'T', 1. )
1393	misfdep(:,:) = INT( zbathy(:,:) )
1394	CALL lbc_lnk( risfdep, 'T', 1. )
1395	CALL lbc_lnk( bathy, 'T', 1. )
1396	zbathy(:,:) = FLOAT( mbathy(:,:) )
1397	CALL lbc_lnk( zbathy, 'T', 1. )
1398	mbathy(:,:) = INT( zbathy(:,:) )
1399	ENDIF
1400	! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1401	DO jj = 1, jpjm1
1402	DO ji = 1, jpim1
1403	IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
1404	IF ( .NOT. ln_iscpl ) THEN
1405	zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
1406	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1407	zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) &
1408	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))
1409	IF (zbathydiff .LE. zrisfdepdiff) THEN
1410	mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1411	bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1412	ELSE
1413	misfdep(ji,jj) = misfdep(ji,jj) - 1
1414	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
1415	END IF
1416	ELSE
1417	misfdep(ji,jj) = misfdep(ji,jj) - 1
1418	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
1419	END IF
1420	ENDIF
1421	END DO
1422	END DO
1423
1424
1425	IF( lk_mpp ) THEN
1426	zbathy(:,:) = FLOAT( misfdep(:,:) )
1427	CALL lbc_lnk( zbathy, 'T', 1. )
1428	misfdep(:,:) = INT( zbathy(:,:) )
1429	CALL lbc_lnk( risfdep, 'T', 1. )
1430	CALL lbc_lnk( bathy, 'T', 1. )
1431	zbathy(:,:) = FLOAT( mbathy(:,:) )
1432	CALL lbc_lnk( zbathy, 'T', 1. )
1433	mbathy(:,:) = INT( zbathy(:,:) )
1434	ENDIF
1435
1436	! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1437	DO jj = 1, jpjm1
1438	DO ji = 1, jpim1
1439	IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1440	IF ( .NOT. ln_iscpl ) THEN
1441	zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
1442	& + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))
1443	zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
1444	& - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1445	IF (zbathydiff .LE. zrisfdepdiff) THEN
1446	mbathy(ji,jj) = mbathy(ji,jj) + 1
1447	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1448	ELSE
1449	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1450	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1451	END IF
1452	ELSE
1453	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1454	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1455	ENDIF
1456	ENDIF
1457	ENDDO
1458	ENDDO
1459
1460	IF( lk_mpp ) THEN
1461	zbathy(:,:) = FLOAT( misfdep(:,:) )
1462	CALL lbc_lnk( zbathy, 'T', 1. )
1463	misfdep(:,:) = INT( zbathy(:,:) )
1464	CALL lbc_lnk( risfdep, 'T', 1. )
1465	CALL lbc_lnk( bathy, 'T', 1. )
1466	zbathy(:,:) = FLOAT( mbathy(:,:) )
1467	CALL lbc_lnk( zbathy, 'T', 1. )
1468	mbathy(:,:) = INT( zbathy(:,:) )
1469	ENDIF
1470
1471	! point U misfdep(ji,jj) EQ bathy(ji,jj+1)
1472	DO jj = 1, jpjm1
1473	DO ji = 1, jpim1
1474	IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1475	IF ( .NOT. ln_iscpl ) THEN
1476	zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
1477	& + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1478	zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) &
1479	& - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))
1480	IF (zbathydiff .LE. zrisfdepdiff) THEN
1481	mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1
1482	bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1483	ELSE
1484	misfdep(ji,jj) = misfdep(ji ,jj) - 1
1485	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
1486	END IF
1487	ELSE
1488	misfdep(ji,jj) = misfdep(ji ,jj) - 1
1489	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
1490	ENDIF
1491	ENDIF
1492	ENDDO
1493	ENDDO
1494
1495	IF( lk_mpp ) THEN
1496	zbathy(:,:) = FLOAT( misfdep(:,:) )
1497	CALL lbc_lnk( zbathy, 'T', 1. )
1498	misfdep(:,:) = INT( zbathy(:,:) )
1499	CALL lbc_lnk( risfdep, 'T', 1. )
1500	CALL lbc_lnk( bathy, 'T', 1. )
1501	zbathy(:,:) = FLOAT( mbathy(:,:) )
1502	CALL lbc_lnk( zbathy, 'T', 1. )
1503	mbathy(:,:) = INT( zbathy(:,:) )
1504	ENDIF
1505	END DO
1506	! end dig bathy/ice shelf to be compatible
1507	! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1508	DO jl = 1,20
1509
1510	! remove single point "bay" on isf coast line in the ice shelf draft'
1511	DO jk = 2, jpk
1512	WHERE (misfdep==0) misfdep=jpk
1513	zmask=0
1514	WHERE (misfdep .LE. jk) zmask=1
1515	DO jj = 2, jpjm1
1516	DO ji = 2, jpim1
1517	IF (misfdep(ji,jj) .EQ. jk) THEN
1518	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1519	IF (ibtest .LE. 1) THEN
1520	risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1521	IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
1522	END IF
1523	END IF
1524	END DO
1525	END DO
1526	END DO
1527	WHERE (misfdep==jpk)
1528	misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1529	END WHERE
1530	IF( lk_mpp ) THEN
1531	zbathy(:,:) = FLOAT( misfdep(:,:) )
1532	CALL lbc_lnk( zbathy, 'T', 1. )
1533	misfdep(:,:) = INT( zbathy(:,:) )
1534	CALL lbc_lnk( risfdep, 'T', 1. )
1535	CALL lbc_lnk( bathy, 'T', 1. )
1536	zbathy(:,:) = FLOAT( mbathy(:,:) )
1537	CALL lbc_lnk( zbathy, 'T', 1. )
1538	mbathy(:,:) = INT( zbathy(:,:) )
1539	ENDIF
1540
1541	! remove single point "bay" on bathy coast line beneath an ice shelf'
1542	DO jk = jpk,1,-1
1543	zmask=0
1544	WHERE (mbathy .GE. jk ) zmask=1
1545	DO jj = 2, jpjm1
1546	DO ji = 2, jpim1
1547	IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
1548	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1549	IF (ibtest .LE. 1) THEN
1550	bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1551	IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
1552	END IF
1553	END IF
1554	END DO
1555	END DO
1556	END DO
1557	WHERE (mbathy==0)
1558	misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1559	END WHERE
1560	IF( lk_mpp ) THEN
1561	zbathy(:,:) = FLOAT( misfdep(:,:) )
1562	CALL lbc_lnk( zbathy, 'T', 1. )
1563	misfdep(:,:) = INT( zbathy(:,:) )
1564	CALL lbc_lnk( risfdep, 'T', 1. )
1565	CALL lbc_lnk( bathy, 'T', 1. )
1566	zbathy(:,:) = FLOAT( mbathy(:,:) )
1567	CALL lbc_lnk( zbathy, 'T', 1. )
1568	mbathy(:,:) = INT( zbathy(:,:) )
1569	ENDIF
1570
1571	! fill hole in ice shelf
1572	zmisfdep = misfdep
1573	zrisfdep = risfdep
1574	WHERE (zmisfdep .LE. 1) zmisfdep=jpk
1575	DO jj = 2, jpjm1
1576	DO ji = 2, jpim1
1577	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
1578	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
1579	IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk
1580	IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk
1581	IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk
1582	IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk
1583	ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1584	IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
1585	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1586	END IF
1587	IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
1588	misfdep(ji,jj) = ibtest
1589	risfdep(ji,jj) = gdepw_1d(ibtest)
1590	ENDIF
1591	ENDDO
1592	ENDDO
1593
1594	IF( lk_mpp ) THEN
1595	zbathy(:,:) = FLOAT( misfdep(:,:) )
1596	CALL lbc_lnk( zbathy, 'T', 1. )
1597	misfdep(:,:) = INT( zbathy(:,:) )
1598	CALL lbc_lnk( risfdep, 'T', 1. )
1599	CALL lbc_lnk( bathy, 'T', 1. )
1600	zbathy(:,:) = FLOAT( mbathy(:,:) )
1601	CALL lbc_lnk( zbathy, 'T', 1. )
1602	mbathy(:,:) = INT( zbathy(:,:) )
1603	ENDIF
1604	!
1605	!! fill hole in bathymetry
1606	zmbathy (:,:)=mbathy (:,:)
1607	DO jj = 2, jpjm1
1608	DO ji = 2, jpim1
1609	ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )
1610	ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)
1611	IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0
1612	IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0
1613	IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0
1614	IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0
1615	ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1616	IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
1617	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1618	END IF
1619	IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
1620	mbathy(ji,jj) = ibtest
1621	bathy(ji,jj) = gdepw_1d(ibtest+1)
1622	ENDIF
1623	END DO
1624	END DO
1625	IF( lk_mpp ) THEN
1626	zbathy(:,:) = FLOAT( misfdep(:,:) )
1627	CALL lbc_lnk( zbathy, 'T', 1. )
1628	misfdep(:,:) = INT( zbathy(:,:) )
1629	CALL lbc_lnk( risfdep, 'T', 1. )
1630	CALL lbc_lnk( bathy, 'T', 1. )
1631	zbathy(:,:) = FLOAT( mbathy(:,:) )
1632	CALL lbc_lnk( zbathy, 'T', 1. )
1633	mbathy(:,:) = INT( zbathy(:,:) )
1634	ENDIF
1635	! if not compatible after all check (ie U point water column less than 2 cells), mask U
1636	DO jj = 1, jpjm1
1637	DO ji = 1, jpim1
1638	IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1639	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
1640	END IF
1641	END DO
1642	END DO
1643	IF( lk_mpp ) THEN
1644	zbathy(:,:) = FLOAT( misfdep(:,:) )
1645	CALL lbc_lnk( zbathy, 'T', 1. )
1646	misfdep(:,:) = INT( zbathy(:,:) )
1647	CALL lbc_lnk( risfdep, 'T', 1. )
1648	CALL lbc_lnk( bathy, 'T', 1. )
1649	zbathy(:,:) = FLOAT( mbathy(:,:) )
1650	CALL lbc_lnk( zbathy, 'T', 1. )
1651	mbathy(:,:) = INT( zbathy(:,:) )
1652	ENDIF
1653	! if not compatible after all check (ie U point water column less than 2 cells), mask U
1654	DO jj = 1, jpjm1
1655	DO ji = 1, jpim1
1656	IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1657	mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;
1658	END IF
1659	END DO
1660	END DO
1661	IF( lk_mpp ) THEN
1662	zbathy(:,:) = FLOAT( misfdep(:,:) )
1663	CALL lbc_lnk( zbathy, 'T', 1. )
1664	misfdep(:,:) = INT( zbathy(:,:) )
1665	CALL lbc_lnk( risfdep, 'T', 1. )
1666	CALL lbc_lnk( bathy, 'T', 1. )
1667	zbathy(:,:) = FLOAT( mbathy(:,:) )
1668	CALL lbc_lnk( zbathy, 'T', 1. )
1669	mbathy(:,:) = INT( zbathy(:,:) )
1670	ENDIF
1671	! if not compatible after all check (ie V point water column less than 2 cells), mask V
1672	DO jj = 1, jpjm1
1673	DO ji = 1, jpi
1674	IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1675	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
1676	END IF
1677	END DO
1678	END DO
1679	IF( lk_mpp ) THEN
1680	zbathy(:,:) = FLOAT( misfdep(:,:) )
1681	CALL lbc_lnk( zbathy, 'T', 1. )
1682	misfdep(:,:) = INT( zbathy(:,:) )
1683	CALL lbc_lnk( risfdep, 'T', 1. )
1684	CALL lbc_lnk( bathy, 'T', 1. )
1685	zbathy(:,:) = FLOAT( mbathy(:,:) )
1686	CALL lbc_lnk( zbathy, 'T', 1. )
1687	mbathy(:,:) = INT( zbathy(:,:) )
1688	ENDIF
1689	! if not compatible after all check (ie V point water column less than 2 cells), mask V
1690	DO jj = 1, jpjm1
1691	DO ji = 1, jpi
1692	IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1693	mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
1694	END IF
1695	END DO
1696	END DO
1697	IF( lk_mpp ) THEN
1698	zbathy(:,:) = FLOAT( misfdep(:,:) )
1699	CALL lbc_lnk( zbathy, 'T', 1. )
1700	misfdep(:,:) = INT( zbathy(:,:) )
1701	CALL lbc_lnk( risfdep, 'T', 1. )
1702	CALL lbc_lnk( bathy, 'T', 1. )
1703	zbathy(:,:) = FLOAT( mbathy(:,:) )
1704	CALL lbc_lnk( zbathy, 'T', 1. )
1705	mbathy(:,:) = INT( zbathy(:,:) )
1706	ENDIF
1707	! if not compatible after all check, mask T
1708	DO jj = 1, jpj
1709	DO ji = 1, jpi
1710	IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1711	misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;
1712	END IF
1713	END DO
1714	END DO
1715
1716	WHERE (mbathy(:,:) == 1)
1717	mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1718	END WHERE
1719	END DO
1720	! end check compatibility ice shelf/bathy
1721	! remove very shallow ice shelf (less than ~ 10m if 75L)
1722	IF( icompt == 0 ) THEN
1723	IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'
1724	ELSE
1725	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'
1726	ENDIF
1727
1728	! compute scale factor and depth at T- and W- points
1729	DO jj = 1, jpj
1730	DO ji = 1, jpi
1731	ik = mbathy(ji,jj)
1732	IF( ik > 0 ) THEN ! ocean point only
1733	! max ocean level case
1734	IF( ik == jpkm1 ) THEN
1735	zdepwp = bathy(ji,jj)
1736	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1737	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1738	e3t_0(ji,jj,ik ) = ze3tp
1739	e3t_0(ji,jj,ik+1) = ze3tp
1740	e3w_0(ji,jj,ik ) = ze3wp
1741	e3w_0(ji,jj,ik+1) = ze3tp
1742	gdepw_0(ji,jj,ik+1) = zdepwp
1743	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1744	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1745	!
1746	ELSE ! standard case
1747	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1748	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1749	ENDIF
1750	! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1751	!gm Bug? check the gdepw_1d
1752	! ... on ik
1753	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1754	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1755	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1756	e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik )
1757	e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1)
1758	! ... on ik+1
1759	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1760	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1761	ENDIF
1762	ENDIF
1763	END DO
1764	END DO
1765	!
1766	it = 0
1767	DO jj = 1, jpj
1768	DO ji = 1, jpi
1769	ik = mbathy(ji,jj)
1770	IF( ik > 0 ) THEN ! ocean point only
1771	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1772	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1773	! test
1774	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1775	IF( zdiff <= 0._wp .AND. lwp ) THEN
1776	it = it + 1
1777	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1778	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1779	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1780	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1781	ENDIF
1782	ENDIF
1783	END DO
1784	END DO
1785	!
1786	! (ISF) Definition of e3t, u, v, w for ISF case
1787	DO jj = 1, jpj
1788	DO ji = 1, jpi
1789	ik = misfdep(ji,jj)
1790	IF( ik > 1 ) THEN ! ice shelf point only
1791	IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
1792	gdepw_0(ji,jj,ik) = risfdep(ji,jj)
1793	!gm Bug? check the gdepw_0
1794	! ... on ik
1795	gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
1796	& * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
1797	& / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
1798	e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
1799	e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1800
1801	IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
1802	e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
1803	ENDIF
1804	! ... on ik / ik-1
1805	e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik)
1806	e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1807	! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1808	gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1809	ENDIF
1810	END DO
1811	END DO
1812
1813	it = 0
1814	DO jj = 1, jpj
1815	DO ji = 1, jpi
1816	ik = misfdep(ji,jj)
1817	IF( ik > 1 ) THEN ! ice shelf point only
1818	e3tp (ji,jj) = e3t_0(ji,jj,ik )
1819	e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
1820	! test
1821	zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
1822	IF( zdiff <= 0. .AND. lwp ) THEN
1823	it = it + 1
1824	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1825	WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
1826	WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1827	WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
1828	ENDIF
1829	ENDIF
1830	END DO
1831	END DO
1832
1833	CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1834	CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1835
1836	IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')
1837
1838	END SUBROUTINE
1839
1840	SUBROUTINE zgr_sco
1841	!!----------------------------------------------------------------------
1842	!! * ROUTINE zgr_sco *
1843	!!
1844	!! ** Purpose : define the s-coordinate system
1845	!!
1846	!! ** Method : s-coordinate
1847	!! The depth of model levels is defined as the product of an
1848	!! analytical function by the local bathymetry, while the vertical
1849	!! scale factors are defined as the product of the first derivative
1850	!! of the analytical function by the bathymetry.
1851	!! (this solution save memory as depth and scale factors are not
1852	!! 3d fields)
1853	!! - Read bathymetry (in meters) at t-point and compute the
1854	!! bathymetry at u-, v-, and f-points.
1855	!! hbatu = mi( hbatt )
1856	!! hbatv = mj( hbatt )
1857	!! hbatf = mi( mj( hbatt ) )
1858	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1859	!! function and its derivative given as function.
1860	!! z_gsigt(k) = fssig (k )
1861	!! z_gsigw(k) = fssig (k-0.5)
1862	!! z_esigt(k) = fsdsig(k )
1863	!! z_esigw(k) = fsdsig(k-0.5)
1864	!! Three options for stretching are give, and they can be modified
1865	!! following the users requirements. Nevertheless, the output as
1866	!! well as the way to compute the model levels and scale factors
1867	!! must be respected in order to insure second order accuracy
1868	!! schemes.
1869	!!
1870	!! The three methods for stretching available are:
1871	!!
1872	!! s_sh94 (Song and Haidvogel 1994)
1873	!! a sinh/tanh function that allows sigma and stretched sigma
1874	!!
1875	!! s_sf12 (Siddorn and Furner 2012?)
1876	!! allows the maintenance of fixed surface and or
1877	!! bottom cell resolutions (cf. geopotential coordinates)
1878	!! within an analytically derived stretched S-coordinate framework.
1879	!!
1880	!! s_tanh (Madec et al 1996)
1881	!! a cosh/tanh function that gives stretched coordinates
1882	!!
1883	!!----------------------------------------------------------------------
1884	!
1885	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1886	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1887	INTEGER :: ios ! Local integer output status for namelist read
1888	REAL(wp) :: zrmax, ztaper ! temporary scalars
1889	REAL(wp) :: zrfact
1890	!
1891	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1892	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1893
1894	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1895	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1896	!!----------------------------------------------------------------------
1897	!
1898	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1899	!
1900	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1901	!
1902	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1903	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1904	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1905
1906	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1907	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1908	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1909	IF(lwm) WRITE ( numond, namzgr_sco )
1910
1911	IF(lwp) THEN ! control print
1912	WRITE(numout,*)
1913	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1914	WRITE(numout,*) '~~~~~~~~~~~'
1915	WRITE(numout,*) ' Namelist namzgr_sco'
1916	WRITE(numout,*) ' stretching coeffs '
1917	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1918	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1919	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1920	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1921	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1922	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1923	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1924	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1925	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1926	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1927	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1928	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1929	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1930	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1931	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1932	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1933	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1934	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1935	ENDIF
1936
1937	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1938	hifu(:,:) = rn_sbot_min
1939	hifv(:,:) = rn_sbot_min
1940	hiff(:,:) = rn_sbot_min
1941
1942	! ! set maximum ocean depth
1943	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1944
1945	DO jj = 1, jpj
1946	DO ji = 1, jpi
1947	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1948	END DO
1949	END DO
1950	! ! =============================
1951	! ! Define the envelop bathymetry (hbatt)
1952	! ! =============================
1953	! use r-value to create hybrid coordinates
1954	zenv(:,:) = bathy(:,:)
1955	!
1956	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1957	DO jj = 1, jpj
1958	DO ji = 1, jpi
1959	IF( bathy(ji,jj) == 0._wp ) THEN
1960	iip1 = MIN( ji+1, jpi )
1961	ijp1 = MIN( jj+1, jpj )
1962	iim1 = MAX( ji-1, 1 )
1963	ijm1 = MAX( jj-1, 1 )
1964	IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + &
1965	& bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1966	zenv(ji,jj) = rn_sbot_min
1967	ENDIF
1968	ENDIF
1969	END DO
1970	END DO
1971	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1972	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1973	!
1974	! smooth the bathymetry (if required)
1975	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1976	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1977	!
1978	jl = 0
1979	zrmax = 1._wp
1980	!
1981	!
1982	! set scaling factor used in reducing vertical gradients
1983	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1984	!
1985	! initialise temporary evelope depth arrays
1986	ztmpi1(:,:) = zenv(:,:)
1987	ztmpi2(:,:) = zenv(:,:)
1988	ztmpj1(:,:) = zenv(:,:)
1989	ztmpj2(:,:) = zenv(:,:)
1990	!
1991	! initialise temporary r-value arrays
1992	zri(:,:) = 1._wp
1993	zrj(:,:) = 1._wp
1994	! ! ================ !
1995	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1996	! ! ================ !
1997	jl = jl + 1
1998	zrmax = 0._wp
1999	! we set zrmax from previous r-values (zri and zrj) first
2000	! if set after current r-value calculation (as previously)
2001	! we could exit DO WHILE prematurely before checking r-value
2002	! of current zenv
2003	DO jj = 1, nlcj
2004	DO ji = 1, nlci
2005	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2006	END DO
2007	END DO
2008	zri(:,:) = 0._wp
2009	zrj(:,:) = 0._wp
2010	DO jj = 1, nlcj
2011	DO ji = 1, nlci
2012	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
2013	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
2014	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2015	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
2016	END IF
2017	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2018	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
2019	END IF
2020	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
2021	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
2022	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
2023	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
2024	END DO
2025	END DO
2026	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
2027	!
2028	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
2029	!
2030	DO jj = 1, nlcj
2031	DO ji = 1, nlci
2032	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2033	END DO
2034	END DO
2035	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
2036	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2037	! ! ================ !
2038	END DO ! End loop !
2039	! ! ================ !
2040	DO jj = 1, jpj
2041	DO ji = 1, jpi
2042	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2043	END DO
2044	END DO
2045	!
2046	! Envelope bathymetry saved in hbatt
2047	hbatt(:,:) = zenv(:,:)
2048	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2049	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2050	DO jj = 1, jpj
2051	DO ji = 1, jpi
2052	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2053	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2054	END DO
2055	END DO
2056	ENDIF
2057	!
2058	IF(lwp) THEN ! Control print
2059	WRITE(numout,*)
2060	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
2061	WRITE(numout,*)
2062	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
2063	IF( nprint == 1 ) THEN
2064	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
2065	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
2066	ENDIF
2067	ENDIF
2068
2069	! ! ==============================
2070	! ! hbatu, hbatv, hbatf fields
2071	! ! ==============================
2072	IF(lwp) THEN
2073	WRITE(numout,*)
2074	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2075	ENDIF
2076	hbatu(:,:) = rn_sbot_min
2077	hbatv(:,:) = rn_sbot_min
2078	hbatf(:,:) = rn_sbot_min
2079	DO jj = 1, jpjm1
2080	DO ji = 1, jpim1 ! NO vector opt.
2081	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
2082	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
2083	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
2084	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2085	END DO
2086	END DO
2087	!
2088	! Apply lateral boundary condition
2089	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2090	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
2091	DO jj = 1, jpj
2092	DO ji = 1, jpi
2093	IF( hbatu(ji,jj) == 0._wp ) THEN
2094	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
2095	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
2096	ENDIF
2097	END DO
2098	END DO
2099	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
2100	DO jj = 1, jpj
2101	DO ji = 1, jpi
2102	IF( hbatv(ji,jj) == 0._wp ) THEN
2103	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
2104	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
2105	ENDIF
2106	END DO
2107	END DO
2108	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
2109	DO jj = 1, jpj
2110	DO ji = 1, jpi
2111	IF( hbatf(ji,jj) == 0._wp ) THEN
2112	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
2113	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
2114	ENDIF
2115	END DO
2116	END DO
2117
2118	!!bug: key_helsinki a verifer
2119	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2120	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2121	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2122	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2123
2124	IF( nprint == 1 .AND. lwp ) THEN
2125	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
2126	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2127	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
2128	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2129	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
2130	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2131	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
2132	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2133	ENDIF
2134	!! helsinki
2135
2136	! ! =======================
2137	! ! s-ccordinate fields (gdep., e3.)
2138	! ! =======================
2139	!
2140	! non-dimensional "sigma" for model level depth at w- and t-levels
2141
2142
2143	!========================================================================
2144	! Song and Haidvogel 1994 (ln_s_sh94=T)
2145	! Siddorn and Furner 2012 (ln_sf12=T)
2146	! or tanh function (both false)
2147	!========================================================================
2148	IF ( ln_s_sh94 ) THEN
2149	CALL s_sh94()
2150	ELSE IF ( ln_s_sf12 ) THEN
2151	CALL s_sf12()
2152	ELSE
2153	CALL s_tanh()
2154	ENDIF
2155
2156	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2157	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2158	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2159	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2160	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2161	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2162	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2163
2164	fsdepw(:,:,:) = gdepw_0 (:,:,:)
2165	fsde3w(:,:,:) = gdep3w_0(:,:,:)
2166	!
2167	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
2168	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
2169	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
2170	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
2171	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
2172	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
2173	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
2174
2175	#if defined key_agrif
2176	! Ensure meaningful vertical scale factors in ghost lines/columns
2177	IF( .NOT. Agrif_Root() ) THEN
2178	!
2179	IF((nbondi == -1).OR.(nbondi == 2)) THEN
2180	e3u_0(1,:,:) = e3u_0(2,:,:)
2181	ENDIF
2182	!
2183	IF((nbondi == 1).OR.(nbondi == 2)) THEN
2184	e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2185	ENDIF
2186	!
2187	IF((nbondj == -1).OR.(nbondj == 2)) THEN
2188	e3v_0(:,1,:) = e3v_0(:,2,:)
2189	ENDIF
2190	!
2191	IF((nbondj == 1).OR.(nbondj == 2)) THEN
2192	e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2193	ENDIF
2194	!
2195	ENDIF
2196	#endif
2197
2198	fsdept(:,:,:) = gdept_0 (:,:,:)
2199	fsdepw(:,:,:) = gdepw_0 (:,:,:)
2200	fsde3w(:,:,:) = gdep3w_0(:,:,:)
2201	fse3t (:,:,:) = e3t_0 (:,:,:)
2202	fse3u (:,:,:) = e3u_0 (:,:,:)
2203	fse3v (:,:,:) = e3v_0 (:,:,:)
2204	fse3f (:,:,:) = e3f_0 (:,:,:)
2205	fse3w (:,:,:) = e3w_0 (:,:,:)
2206	fse3uw(:,:,:) = e3uw_0 (:,:,:)
2207	fse3vw(:,:,:) = e3vw_0 (:,:,:)
2208	!!
2209	! HYBRID :
2210	DO jj = 1, jpj
2211	DO ji = 1, jpi
2212	DO jk = 1, jpkm1
2213	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
2214	END DO
2215	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
2216	END DO
2217	END DO
2218	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
2219	& ' MAX ', MAXVAL( mbathy(:,:) )
2220
2221	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
2222	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
2223	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
2224	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
2225	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
2226	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
2227	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
2228	& ' w ', MINVAL( e3w_0 (:,:,:) )
2229
2230	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
2231	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
2232	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
2233	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
2234	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
2235	& ' w ', MAXVAL( e3w_0 (:,:,:) )
2236	ENDIF
2237	! END DO
2238	IF(lwp) THEN ! selected vertical profiles
2239	WRITE(numout,*)
2240	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2241	WRITE(numout,*) ' ~~~~~~ --------------------'
2242	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2243	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
2244	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2245	DO jj = mj0(20), mj1(20)
2246	DO ji = mi0(20), mi1(20)
2247	WRITE(numout,*)
2248	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2249	WRITE(numout,*) ' ~~~~~~ --------------------'
2250	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2251	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2252	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2253	END DO
2254	END DO
2255	DO jj = mj0(74), mj1(74)
2256	DO ji = mi0(100), mi1(100)
2257	WRITE(numout,*)
2258	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2259	WRITE(numout,*) ' ~~~~~~ --------------------'
2260	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2261	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2262	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2263	END DO
2264	END DO
2265	ENDIF
2266
2267	!================================================================================
2268	! check the coordinate makes sense
2269	!================================================================================
2270	DO ji = 1, jpi
2271	DO jj = 1, jpj
2272
2273	IF( hbatt(ji,jj) > 0._wp) THEN
2274	DO jk = 1, mbathy(ji,jj)
2275	! check coordinate is monotonically increasing
2276	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2277	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2278	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2279	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2280	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2281	CALL ctl_stop( ctmp1 )
2282	ENDIF
2283	! and check it has never gone negative
2284	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2285	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2286	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2287	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2288	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2289	CALL ctl_stop( ctmp1 )
2290	ENDIF
2291	! and check it never exceeds the total depth
2292	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2293	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2294	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2295	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2296	CALL ctl_stop( ctmp1 )
2297	ENDIF
2298	END DO
2299
2300	DO jk = 1, mbathy(ji,jj)-1
2301	! and check it never exceeds the total depth
2302	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2303	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2304	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2305	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2306	CALL ctl_stop( ctmp1 )
2307	ENDIF
2308	END DO
2309
2310	ENDIF
2311
2312	END DO
2313	END DO
2314	!
2315	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2316	!
2317	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
2318	!
2319	END SUBROUTINE zgr_sco
2320
2321	!!======================================================================
2322	SUBROUTINE s_sh94()
2323
2324	!!----------------------------------------------------------------------
2325	!! * ROUTINE s_sh94 *
2326	!!
2327	!! ** Purpose : stretch the s-coordinate system
2328	!!
2329	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
2330	!! mixed S/sigma coordinate
2331	!!
2332	!! Reference : Song and Haidvogel 1994.
2333	!!----------------------------------------------------------------------
2334	!
2335	INTEGER :: ji, jj, jk ! dummy loop argument
2336	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2337	!
2338	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2339	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2340
2341	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2342	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2343
2344	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2345	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2346	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2347	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2348
2349	DO ji = 1, jpi
2350	DO jj = 1, jpj
2351
2352	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2353	DO jk = 1, jpk
2354	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2355	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2356	END DO
2357	ELSE ! shallow water, uniform sigma
2358	DO jk = 1, jpk
2359	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
2360	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2361	END DO
2362	ENDIF
2363	!
2364	DO jk = 1, jpkm1
2365	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2366	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2367	END DO
2368	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
2369	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2370	!
2371	! Coefficients for vertical depth as the sum of e3w scale factors
2372	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2373	DO jk = 2, jpk
2374	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2375	END DO
2376	!
2377	DO jk = 1, jpk
2378	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2379	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2380	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
2381	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
2382	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
2383	END DO
2384	!
2385	END DO ! for all jj's
2386	END DO ! for all ji's
2387
2388	DO ji = 1, jpim1
2389	DO jj = 1, jpjm1
2390	DO jk = 1, jpk
2391	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
2392	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2393	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
2394	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2395	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
2396	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
2397	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2398	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
2399	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2400	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
2401	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2402	!
2403	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2404	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2405	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2406	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2407	!
2408	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2409	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2410	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2411	END DO
2412	END DO
2413	END DO
2414
2415	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2416	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2417
2418	END SUBROUTINE s_sh94
2419
2420	SUBROUTINE s_sf12
2421
2422	!!----------------------------------------------------------------------
2423	!! * ROUTINE s_sf12 *
2424	!!
2425	!! ** Purpose : stretch the s-coordinate system
2426	!!
2427	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
2428	!! mixed S/sigma/Z coordinate
2429	!!
2430	!! This method allows the maintenance of fixed surface and or
2431	!! bottom cell resolutions (cf. geopotential coordinates)
2432	!! within an analytically derived stretched S-coordinate framework.
2433	!!
2434	!!
2435	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2436	!!----------------------------------------------------------------------
2437	!
2438	INTEGER :: ji, jj, jk ! dummy loop argument
2439	REAL(wp) :: zsmth ! smoothing around critical depth
2440	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
2441	!
2442	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2443	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2444
2445	!
2446	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2447	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2448
2449	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2450	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2451	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2452	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2453
2454	DO ji = 1, jpi
2455	DO jj = 1, jpj
2456
2457	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2458
2459	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
2460	! could be changed by users but care must be taken to do so carefully
2461	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2462
2463	zzs = rn_zs / hbatt(ji,jj)
2464
2465	IF (rn_efold /= 0.0_wp) THEN
2466	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2467	ELSE
2468	zsmth = 1.0_wp
2469	ENDIF
2470
2471	DO jk = 1, jpk
2472	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2473	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2474	ENDDO
2475	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
2476	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
2477
2478	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2479
2480	DO jk = 1, jpk
2481	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2482	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2483	END DO
2484
2485	ELSE ! shallow water, z coordinates
2486
2487	DO jk = 1, jpk
2488	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2489	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2490	END DO
2491
2492	ENDIF
2493
2494	DO jk = 1, jpkm1
2495	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2496	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2497	END DO
2498	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
2499	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2500
2501	! Coefficients for vertical depth as the sum of e3w scale factors
2502	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2503	DO jk = 2, jpk
2504	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2505	END DO
2506
2507	DO jk = 1, jpk
2508	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2509	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2510	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2511	END DO
2512
2513	ENDDO ! for all jj's
2514	ENDDO ! for all ji's
2515
2516	DO ji=1,jpi-1
2517	DO jj=1,jpj-1
2518
2519	DO jk = 1, jpk
2520	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
2521	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2522	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
2523	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2524	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
2525	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
2526	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2527	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
2528	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2529	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
2530	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2531
2532	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2533	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2534	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2535	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2536	!
2537	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2538	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2539	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2540	END DO
2541
2542	ENDDO
2543	ENDDO
2544	!
2545	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2546	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2547	CALL lbc_lnk(e3w_0 ,'T',1.)
2548	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2549	!
2550	! ! =============
2551
2552	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2553	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2554
2555	END SUBROUTINE s_sf12
2556
2557	SUBROUTINE s_tanh()
2558
2559	!!----------------------------------------------------------------------
2560	!! * ROUTINE s_tanh*
2561	!!
2562	!! ** Purpose : stretch the s-coordinate system
2563	!!
2564	!! ** Method : s-coordinate stretch
2565	!!
2566	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2567	!!----------------------------------------------------------------------
2568
2569	INTEGER :: ji, jj, jk ! dummy loop argument
2570	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2571
2572	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2573	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2574
2575	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2576	CALL wrk_alloc( jpk, z_esigt, z_esigw )
2577
2578	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
2579	z_esigt = 0._wp ; z_esigw = 0._wp
2580
2581	DO jk = 1, jpk
2582	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2583	z_gsigt(jk) = -fssig( REAL(jk,wp) )
2584	END DO
2585	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
2586	!
2587	! Coefficients for vertical scale factors at w-, t- levels
2588	!!gm bug : define it from analytical function, not like juste bellow....
2589	!!gm or betteroffer the 2 possibilities....
2590	DO jk = 1, jpkm1
2591	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
2592	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2593	END DO
2594	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
2595	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2596	!
2597	! Coefficients for vertical depth as the sum of e3w scale factors
2598	z_gsi3w(1) = 0.5_wp * z_esigw(1)
2599	DO jk = 2, jpk
2600	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2601	END DO
2602	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2603	DO jk = 1, jpk
2604	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2605	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2606	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
2607	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
2608	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
2609	END DO
2610	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2611	DO jj = 1, jpj
2612	DO ji = 1, jpi
2613	DO jk = 1, jpk
2614	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2615	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2616	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2617	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2618	!
2619	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2620	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2621	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2622	END DO
2623	END DO
2624	END DO
2625
2626	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2627	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2628
2629	END SUBROUTINE s_tanh
2630
2631	FUNCTION fssig( pk ) RESULT( pf )
2632	!!----------------------------------------------------------------------
2633	!! * ROUTINE fssig *
2634	!!
2635	!! ** Purpose : provide the analytical function in s-coordinate
2636	!!
2637	!! ** Method : the function provide the non-dimensional position of
2638	!! T and W (i.e. between 0 and 1)
2639	!! T-points at integer values (between 1 and jpk)
2640	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2641	!!----------------------------------------------------------------------
2642	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2643	REAL(wp) :: pf ! sigma value
2644	!!----------------------------------------------------------------------
2645	!
2646	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2647	& - TANH( rn_thetb * rn_theta ) ) &
2648	& * ( COSH( rn_theta ) &
2649	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2650	& / ( 2._wp * SINH( rn_theta ) )
2651	!
2652	END FUNCTION fssig
2653
2654
2655	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2656	!!----------------------------------------------------------------------
2657	!! * ROUTINE fssig1 *
2658	!!
2659	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2660	!!
2661	!! ** Method : the function provides the non-dimensional position of
2662	!! T and W (i.e. between 0 and 1)
2663	!! T-points at integer values (between 1 and jpk)
2664	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2665	!!----------------------------------------------------------------------
2666	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2667	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2668	REAL(wp) :: pf1 ! sigma value
2669	!!----------------------------------------------------------------------
2670	!
2671	IF ( rn_theta == 0 ) then ! uniform sigma
2672	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2673	ELSE ! stretched sigma
2674	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2675	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2676	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2677	ENDIF
2678	!
2679	END FUNCTION fssig1
2680
2681
2682	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2683	!!----------------------------------------------------------------------
2684	!! * ROUTINE fgamma *
2685	!!
2686	!! ** Purpose : provide analytical function for the s-coordinate
2687	!!
2688	!! ** Method : the function provides the non-dimensional position of
2689	!! T and W (i.e. between 0 and 1)
2690	!! T-points at integer values (between 1 and jpk)
2691	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2692	!!
2693	!! This method allows the maintenance of fixed surface and or
2694	!! bottom cell resolutions (cf. geopotential coordinates)
2695	!! within an analytically derived stretched S-coordinate framework.
2696	!!
2697	!! Reference : Siddorn and Furner, in prep
2698	!!----------------------------------------------------------------------
2699	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2700	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2701	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2702	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2703	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2704	REAL(wp) :: za1,za2,za3 ! local variables
2705	REAL(wp) :: zn1,zn2 ! local variables
2706	REAL(wp) :: za,zb,zx ! local variables
2707	integer :: jk
2708	!!----------------------------------------------------------------------
2709	!
2710
2711	zn1 = 1./(jpk-1.)
2712	zn2 = 1. - zn1
2713
2714	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2715	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2716	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2717
2718	za = pzb - za3*(pzs-za1)-za2
2719	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2720	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2721	zx = 1.0_wp-za/2.0_wp-zb
2722
2723	DO jk = 1, jpk
2724	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2725	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2726	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2727	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2728	ENDDO
2729
2730	!
2731	END FUNCTION fgamma
2732
2733	!!======================================================================
2734	END MODULE domzgr

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

Download in other formats: