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 @ 5920

Last change on this file since 5920 was 5920, checked in by mathiot, 8 years ago
ice sheet coupling: add treshold to define grounded area, remove useless va riable, change some variable name + add some namelist parameter
Property svn:keywords set to `Id`
File size: 135.1 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
532	! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
533	WHERE ( bathy(:,:) .LE. risfdep(:,:) + rn_isfhmin )
534	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
535	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
536	END WHERE
537	END IF
538	!
539	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
540	!
541	IF( nn_cla == 0 ) THEN
542	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
543	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
544	DO ji = mi0(ii0), mi1(ii1)
545	DO jj = mj0(ij0), mj1(ij1)
546	bathy(ji,jj) = 284._wp
547	END DO
548	END DO
549	IF(lwp) WRITE(numout,*)
550	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
551	!
552	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
553	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
554	DO ji = mi0(ii0), mi1(ii1)
555	DO jj = mj0(ij0), mj1(ij1)
556	bathy(ji,jj) = 137._wp
557	END DO
558	END DO
559	IF(lwp) WRITE(numout,*)
560	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
561	ENDIF
562	!
563	ENDIF
564	!
565	ENDIF
566	! ! =============== !
567	ELSE ! error !
568	! ! =============== !
569	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
570	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
571	ENDIF
572	!
573	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
574	!
575	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
576	! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
577	IF ( ln_isfcav ) THEN
578	WHERE ( bathy == risfdep )
579	bathy = 0.0_wp ; risfdep = 0.0_wp
580	END WHERE
581	END IF
582	! end patch
583	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
584	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
585	ENDIF
586	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
587	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
588	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
589	END WHERE
590	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
591	ENDIF
592	!
593	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')
594	!
595	END SUBROUTINE zgr_bat
596
597
598	SUBROUTINE zgr_bat_zoom
599	!!----------------------------------------------------------------------
600	!! * ROUTINE zgr_bat_zoom *
601	!!
602	!! ** Purpose : - Close zoom domain boundary if necessary
603	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
604	!!
605	!! ** Method :
606	!!
607	!! ** Action : - update mbathy: level bathymetry (in level index)
608	!!----------------------------------------------------------------------
609	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
610	!!----------------------------------------------------------------------
611	!
612	IF(lwp) WRITE(numout,*)
613	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
614	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
615	!
616	! Zoom domain
617	! ===========
618	!
619	! Forced closed boundary if required
620	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
621	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
622	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
623	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
624	!
625	! Configuration specific domain modifications
626	! (here, ORCA arctic configuration: suppress Med Sea)
627	IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
628	SELECT CASE ( jp_cfg )
629	! ! =======================
630	CASE ( 2 ) ! ORCA_R2 configuration
631	! ! =======================
632	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
633	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
634	ij0 = 98 ; ij1 = 110
635	! ! =======================
636	CASE ( 05 ) ! ORCA_R05 configuration
637	! ! =======================
638	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
639	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
640	ij0 = 314 ; ij1 = 370
641	END SELECT
642	!
643	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
644	!
645	ENDIF
646	!
647	END SUBROUTINE zgr_bat_zoom
648
649
650	SUBROUTINE zgr_bat_ctl
651	!!----------------------------------------------------------------------
652	!! * ROUTINE zgr_bat_ctl *
653	!!
654	!! ** Purpose : check the bathymetry in levels
655	!!
656	!! ** Method : The array mbathy is checked to verified its consistency
657	!! with the model options. in particular:
658	!! mbathy must have at least 1 land grid-points (mbathy<=0)
659	!! along closed boundary.
660	!! mbathy must be cyclic IF jperio=1.
661	!! mbathy must be lower or equal to jpk-1.
662	!! isolated ocean grid points are suppressed from mbathy
663	!! since they are only connected to remaining
664	!! ocean through vertical diffusion.
665	!! C A U T I O N : mbathy will be modified during the initializa-
666	!! tion phase to become the number of non-zero w-levels of a water
667	!! column, with a minimum value of 1.
668	!!
669	!! ** Action : - update mbathy: level bathymetry (in level index)
670	!! - update bathy : meter bathymetry (in meters)
671	!!----------------------------------------------------------------------
672	!!
673	INTEGER :: ji, jj, jl ! dummy loop indices
674	INTEGER :: icompt, ibtest, ikmax ! temporary integers
675	REAL(wp), POINTER, DIMENSION(:,:) :: zbathy
676
677	!!----------------------------------------------------------------------
678	!
679	IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl')
680	!
681	CALL wrk_alloc( jpi, jpj, zbathy )
682	!
683	IF(lwp) WRITE(numout,*)
684	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
685	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
686	! ! Suppress isolated ocean grid points
687	IF(lwp) WRITE(numout,*)
688	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
689	IF(lwp) WRITE(numout,*)' -----------------------------------'
690	icompt = 0
691	DO jl = 1, 2
692	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
693	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
694	mbathy(jpi,:) = mbathy( 2 ,:)
695	ENDIF
696	DO jj = 2, jpjm1
697	DO ji = 2, jpim1
698	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
699	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
700	IF( ibtest < mbathy(ji,jj) ) THEN
701	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
702	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
703	mbathy(ji,jj) = ibtest
704	icompt = icompt + 1
705	ENDIF
706	END DO
707	END DO
708	END DO
709	IF( lk_mpp ) CALL mpp_sum( icompt )
710	IF( icompt == 0 ) THEN
711	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
712	ELSE
713	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
714	ENDIF
715	IF( lk_mpp ) THEN
716	zbathy(:,:) = FLOAT( mbathy(:,:) )
717	CALL lbc_lnk( zbathy, 'T', 1._wp )
718	mbathy(:,:) = INT( zbathy(:,:) )
719	ENDIF
720	! ! East-west cyclic boundary conditions
721	IF( nperio == 0 ) THEN
722	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
723	IF( lk_mpp ) THEN
724	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
725	IF( jperio /= 1 ) mbathy(1,:) = 0
726	ENDIF
727	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
728	IF( jperio /= 1 ) mbathy(nlci,:) = 0
729	ENDIF
730	ELSE
731	IF( ln_zco .OR. ln_zps ) THEN
732	mbathy( 1 ,:) = 0
733	mbathy(jpi,:) = 0
734	ELSE
735	mbathy( 1 ,:) = jpkm1
736	mbathy(jpi,:) = jpkm1
737	ENDIF
738	ENDIF
739	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
740	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
741	mbathy( 1 ,:) = mbathy(jpim1,:)
742	mbathy(jpi,:) = mbathy( 2 ,:)
743	ELSEIF( nperio == 2 ) THEN
744	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
745	ELSE
746	IF(lwp) WRITE(numout,*) ' e r r o r'
747	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
748	! STOP 'dom_mba'
749	ENDIF
750	! Boundary condition on mbathy
751	IF( .NOT.lk_mpp ) THEN
752	!!gm !!bug ??? think about it !
753	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
754	zbathy(:,:) = FLOAT( mbathy(:,:) )
755	CALL lbc_lnk( zbathy, 'T', 1._wp )
756	mbathy(:,:) = INT( zbathy(:,:) )
757	ENDIF
758	! Number of ocean level inferior or equal to jpkm1
759	ikmax = 0
760	DO jj = 1, jpj
761	DO ji = 1, jpi
762	ikmax = MAX( ikmax, mbathy(ji,jj) )
763	END DO
764	END DO
765	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
766	IF( ikmax > jpkm1 ) THEN
767	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
768	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
769	ELSE IF( ikmax < jpkm1 ) THEN
770	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
771	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
772	ENDIF
773
774	IF( lwp .AND. nprint == 1 ) THEN ! control print
775	WRITE(numout,*)
776	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
777	WRITE(numout,*) ' ------------------'
778	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
779	WRITE(numout,*)
780	ENDIF
781	!
782	CALL wrk_dealloc( jpi, jpj, zbathy )
783	!
784	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl')
785	!
786	END SUBROUTINE zgr_bat_ctl
787
788
789	SUBROUTINE zgr_bot_level
790	!!----------------------------------------------------------------------
791	!! * ROUTINE zgr_bot_level *
792	!!
793	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
794	!!
795	!! ** Method : computes from mbathy with a minimum value of 1 over land
796	!!
797	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
798	!! ocean level at t-, u- & v-points
799	!! (min value = 1 over land)
800	!!----------------------------------------------------------------------
801	!!
802	INTEGER :: ji, jj ! dummy loop indices
803	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
804	!!----------------------------------------------------------------------
805	!
806	IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level')
807	!
808	CALL wrk_alloc( jpi, jpj, zmbk )
809	!
810	IF(lwp) WRITE(numout,*)
811	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
812	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
813	!
814	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
815
816	! ! bottom k-index of W-level = mbkt+1
817	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
818	DO ji = 1, jpim1
819	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
820	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
821	END DO
822	END DO
823	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
824	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
825	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
826	!
827	CALL wrk_dealloc( jpi, jpj, zmbk )
828	!
829	IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level')
830	!
831	END SUBROUTINE zgr_bot_level
832
833	SUBROUTINE zgr_top_level
834	!!----------------------------------------------------------------------
835	!! * ROUTINE zgr_bot_level *
836	!!
837	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
838	!!
839	!! ** Method : computes from misfdep with a minimum value of 1
840	!!
841	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
842	!! ocean level at t-, u- & v-points
843	!! (min value = 1)
844	!!----------------------------------------------------------------------
845	!!
846	INTEGER :: ji, jj ! dummy loop indices
847	REAL(wp), POINTER, DIMENSION(:,:) :: zmik
848	!!----------------------------------------------------------------------
849	!
850	IF( nn_timing == 1 ) CALL timing_start('zgr_top_level')
851	!
852	CALL wrk_alloc( jpi, jpj, zmik )
853	!
854	IF(lwp) WRITE(numout,*)
855	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
856	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
857	!
858	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
859	! ! top k-index of W-level (=mikt)
860	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
861	DO ji = 1, jpim1
862	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
863	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
864	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
865	END DO
866	END DO
867
868	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
869	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
870	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
871	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
872	!
873	CALL wrk_dealloc( jpi, jpj, zmik )
874	!
875	IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level')
876	!
877	END SUBROUTINE zgr_top_level
878
879	SUBROUTINE zgr_zco
880	!!----------------------------------------------------------------------
881	!! * ROUTINE zgr_zco *
882	!!
883	!! ** Purpose : define the z-coordinate system
884	!!
885	!! ** Method : set 3D coord. arrays to reference 1D array
886	!!----------------------------------------------------------------------
887	INTEGER :: jk
888	!!----------------------------------------------------------------------
889	!
890	IF( nn_timing == 1 ) CALL timing_start('zgr_zco')
891	!
892	DO jk = 1, jpk
893	gdept_0 (:,:,jk) = gdept_1d(jk)
894	gdepw_0 (:,:,jk) = gdepw_1d(jk)
895	gdep3w_0(:,:,jk) = gdepw_1d(jk)
896	e3t_0 (:,:,jk) = e3t_1d (jk)
897	e3u_0 (:,:,jk) = e3t_1d (jk)
898	e3v_0 (:,:,jk) = e3t_1d (jk)
899	e3f_0 (:,:,jk) = e3t_1d (jk)
900	e3w_0 (:,:,jk) = e3w_1d (jk)
901	e3uw_0 (:,:,jk) = e3w_1d (jk)
902	e3vw_0 (:,:,jk) = e3w_1d (jk)
903	END DO
904	!
905	IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')
906	!
907	END SUBROUTINE zgr_zco
908
909
910	SUBROUTINE zgr_zps
911	!!----------------------------------------------------------------------
912	!! * ROUTINE zgr_zps *
913	!!
914	!! ** Purpose : the depth and vertical scale factor in partial step
915	!! z-coordinate case
916	!!
917	!! ** Method : Partial steps : computes the 3D vertical scale factors
918	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
919	!! a partial step representation of bottom topography.
920	!!
921	!! The reference depth of model levels is defined from an analytical
922	!! function the derivative of which gives the reference vertical
923	!! scale factors.
924	!! From depth and scale factors reference, we compute there new value
925	!! with partial steps on 3d arrays ( i, j, k ).
926	!!
927	!! w-level: gdepw_0(i,j,k) = gdep(k)
928	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
929	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
930	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
931	!!
932	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
933	!! we find the mbathy index of the depth at each grid point.
934	!! This leads us to three cases:
935	!!
936	!! - bathy = 0 => mbathy = 0
937	!! - 1 < mbathy < jpkm1
938	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
939	!!
940	!! Then, for each case, we find the new depth at t- and w- levels
941	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
942	!! and f-points.
943	!!
944	!! This routine is given as an example, it must be modified
945	!! following the user s desiderata. nevertheless, the output as
946	!! well as the way to compute the model levels and scale factors
947	!! must be respected in order to insure second order accuracy
948	!! schemes.
949	!!
950	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
951	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
952	!!
953	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
954	!!----------------------------------------------------------------------
955	!!
956	INTEGER :: ji, jj, jk ! dummy loop indices
957	INTEGER :: ik, it, ikb, ikt ! temporary integers
958	LOGICAL :: ll_print ! Allow control print for debugging
959	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
960	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
961	REAL(wp) :: zdiff ! temporary scalar
962	REAL(wp) :: zmax ! temporary scalar
963	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
964	!!---------------------------------------------------------------------
965	!
966	IF( nn_timing == 1 ) CALL timing_start('zgr_zps')
967	!
968	CALL wrk_alloc( jpi, jpj, jpk, zprt )
969	!
970	IF(lwp) WRITE(numout,*)
971	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
972	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
973	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
974
975	ll_print = .FALSE. ! Local variable for debugging
976
977	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
978	WRITE(numout,*)
979	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
980	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
981	ENDIF
982
983
984	! bathymetry in level (from bathy_meter)
985	! ===================
986	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
987	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
988	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
989	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
990	END WHERE
991
992	! Compute mbathy for ocean points (i.e. the number of ocean levels)
993	! find the number of ocean levels such that the last level thickness
994	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
995	! e3t_1d is the reference level thickness
996	DO jk = jpkm1, 1, -1
997	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
998	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
999	END DO
1000
1001	! Scale factors and depth at T- and W-points
1002	DO jk = 1, jpk ! intitialization to the reference z-coordinate
1003	gdept_0(:,:,jk) = gdept_1d(jk)
1004	gdepw_0(:,:,jk) = gdepw_1d(jk)
1005	e3t_0 (:,:,jk) = e3t_1d (jk)
1006	e3w_0 (:,:,jk) = e3w_1d (jk)
1007	END DO
1008
1009	! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
1010	IF ( ln_isfcav ) CALL zgr_isf
1011
1012	! Scale factors and depth at T- and W-points
1013	IF ( .NOT. ln_isfcav ) THEN
1014	DO jj = 1, jpj
1015	DO ji = 1, jpi
1016	ik = mbathy(ji,jj)
1017	IF( ik > 0 ) THEN ! ocean point only
1018	! max ocean level case
1019	IF( ik == jpkm1 ) THEN
1020	zdepwp = bathy(ji,jj)
1021	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1022	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1023	e3t_0(ji,jj,ik ) = ze3tp
1024	e3t_0(ji,jj,ik+1) = ze3tp
1025	e3w_0(ji,jj,ik ) = ze3wp
1026	e3w_0(ji,jj,ik+1) = ze3tp
1027	gdepw_0(ji,jj,ik+1) = zdepwp
1028	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1029	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1030	!
1031	ELSE ! standard case
1032	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1033	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1034	ENDIF
1035	!gm Bug? check the gdepw_1d
1036	! ... on ik
1037	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1038	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1039	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1040	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1041	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1042	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1043	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1044	! ... on ik+1
1045	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1046	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1047	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1048	ENDIF
1049	ENDIF
1050	END DO
1051	END DO
1052	!
1053	it = 0
1054	DO jj = 1, jpj
1055	DO ji = 1, jpi
1056	ik = mbathy(ji,jj)
1057	IF( ik > 0 ) THEN ! ocean point only
1058	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1059	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1060	! test
1061	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1062	IF( zdiff <= 0._wp .AND. lwp ) THEN
1063	it = it + 1
1064	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1065	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1066	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1067	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1068	ENDIF
1069	ENDIF
1070	END DO
1071	END DO
1072	END IF
1073	!
1074	! Scale factors and depth at U-, V-, UW and VW-points
1075	DO jk = 1, jpk ! initialisation to z-scale factors
1076	e3u_0 (:,:,jk) = e3t_1d(jk)
1077	e3v_0 (:,:,jk) = e3t_1d(jk)
1078	e3uw_0(:,:,jk) = e3w_1d(jk)
1079	e3vw_0(:,:,jk) = e3w_1d(jk)
1080	END DO
1081
1082	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1083	DO jj = 1, jpjm1
1084	DO ji = 1, fs_jpim1 ! vector opt.
1085	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1086	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1087	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1088	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1089	END DO
1090	END DO
1091	END DO
1092	IF ( ln_isfcav ) THEN
1093	! (ISF) define e3uw (adapted for 2 cells in the water column)
1094	DO jj = 2, jpjm1
1095	DO ji = 2, fs_jpim1 ! vector opt.
1096	ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1097	ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1098	IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &
1099	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )
1100	ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1101	ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1102	IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &
1103	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )
1104	END DO
1105	END DO
1106	END IF
1107
1108	CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1109	CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1110	!
1111
1112	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1113	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1114	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1115	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1116	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1117	END DO
1118
1119	! Scale factor at F-point
1120	DO jk = 1, jpk ! initialisation to z-scale factors
1121	e3f_0(:,:,jk) = e3t_1d(jk)
1122	END DO
1123	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1124	DO jj = 1, jpjm1
1125	DO ji = 1, fs_jpim1 ! vector opt.
1126	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1127	END DO
1128	END DO
1129	END DO
1130	CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1131	!
1132	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1133	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1134	END DO
1135	!!gm bug ? : must be a do loop with mj0,mj1
1136	!
1137	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1138	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1139	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1140	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1141	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1142
1143	! Control of the sign
1144	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1145	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1146	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1147	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1148
1149	! Compute gdep3w_0 (vertical sum of e3w)
1150	IF ( ln_isfcav ) THEN ! if cavity
1151	WHERE (misfdep == 0) misfdep = 1
1152	DO jj = 1,jpj
1153	DO ji = 1,jpi
1154	gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1155	DO jk = 2, misfdep(ji,jj)
1156	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1157	END DO
1158	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))
1159	DO jk = misfdep(ji,jj) + 1, jpk
1160	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1161	END DO
1162	END DO
1163	END DO
1164	ELSE ! no cavity
1165	gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1166	DO jk = 2, jpk
1167	gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1168	END DO
1169	END IF
1170	! ! ================= !
1171	IF(lwp .AND. ll_print) THEN ! Control print !
1172	! ! ================= !
1173	DO jj = 1,jpj
1174	DO ji = 1, jpi
1175	ik = MAX( mbathy(ji,jj), 1 )
1176	zprt(ji,jj,1) = e3t_0 (ji,jj,ik)
1177	zprt(ji,jj,2) = e3w_0 (ji,jj,ik)
1178	zprt(ji,jj,3) = e3u_0 (ji,jj,ik)
1179	zprt(ji,jj,4) = e3v_0 (ji,jj,ik)
1180	zprt(ji,jj,5) = e3f_0 (ji,jj,ik)
1181	zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1182	END DO
1183	END DO
1184	WRITE(numout,*)
1185	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1186	WRITE(numout,*)
1187	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1188	WRITE(numout,*)
1189	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1190	WRITE(numout,*)
1191	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1192	WRITE(numout,*)
1193	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1194	WRITE(numout,*)
1195	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1196	ENDIF
1197	!
1198	CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1199	!
1200	IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
1201	!
1202	END SUBROUTINE zgr_zps
1203
1204	SUBROUTINE zgr_isf
1205	!!----------------------------------------------------------------------
1206	!! * ROUTINE zgr_isf *
1207	!!
1208	!! ** Purpose : check the bathymetry in levels
1209	!!
1210	!! ** Method : THe water column have to contained at least 2 cells
1211	!! Bathymetry and isfdraft are modified (dig/close) to respect
1212	!! this criterion.
1213	!!
1214	!!
1215	!! ** Action : - test compatibility between isfdraft and bathy
1216	!! - bathy and isfdraft are modified
1217	!!----------------------------------------------------------------------
1218	!!
1219	INTEGER :: ji, jj, jl, jk ! dummy loop indices
1220	INTEGER :: ik, it ! temporary integers
1221	INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
1222	REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t
1223	REAL(wp) :: zmax ! Maximum and minimum depth
1224	REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar
1225	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1226	REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t
1227	REAL(wp) :: zdiff ! temporary scalar
1228	REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)
1229	INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)
1230	!!---------------------------------------------------------------------
1231	!
1232	IF( nn_timing == 1 ) CALL timing_start('zgr_isf')
1233	!
1234	CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
1235	CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1236
1237
1238	! (ISF) compute misfdep
1239	WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ; misfdep(:,:) = 1 ! open water : set misfdep to 1
1240	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
1241	END WHERE
1242
1243	! Compute misfdep for ocean points (i.e. first wet level)
1244	! find the first ocean level such that the first level thickness
1245	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1246	! e3t_0 is the reference level thickness
1247	DO jk = 2, jpkm1
1248	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1249	WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
1250	END DO
1251	WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
1252	risfdep(:,:) = 0. ; misfdep(:,:) = 1
1253	END WHERE
1254
1255	! 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
1256	icompt = 0
1257	! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
1258	DO jl = 1, 10
1259	! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
1260	WHERE (bathy(:,:) .LE. risfdep(:,:) + rn_isfhmin)
1261	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
1262	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
1263	END WHERE
1264	WHERE (mbathy(:,:) <= 0)
1265	misfdep(:,:) = 0; risfdep(:,:) = 0._wp
1266	mbathy (:,:) = 0; bathy (:,:) = 0._wp
1267	ENDWHERE
1268	IF( lk_mpp ) THEN
1269	zbathy(:,:) = FLOAT( misfdep(:,:) )
1270	CALL lbc_lnk( zbathy, 'T', 1. )
1271	misfdep(:,:) = INT( zbathy(:,:) )
1272	CALL lbc_lnk( risfdep, 'T', 1. )
1273	CALL lbc_lnk( bathy, 'T', 1. )
1274	zbathy(:,:) = FLOAT( mbathy(:,:) )
1275	CALL lbc_lnk( zbathy, 'T', 1. )
1276	mbathy(:,:) = INT( zbathy(:,:) )
1277	ENDIF
1278	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1279	misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west
1280	misfdep(jpi,:) = misfdep( 2 ,:)
1281	ENDIF
1282
1283	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1284	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
1285	mbathy(jpi,:) = mbathy( 2 ,:)
1286	ENDIF
1287
1288	! split last cell if possible (only where water column is 2 cell or less)
1289	! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
1290	IF ( .NOT. ln_iscpl) THEN
1291	DO jk = jpkm1, 1, -1
1292	zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1293	WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
1294	mbathy(:,:) = jk
1295	bathy(:,:) = zmax
1296	END WHERE
1297	END DO
1298	END IF
1299
1300	! split top cell if possible (only where water column is 2 cell or less)
1301	DO jk = 2, jpkm1
1302	zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1303	WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
1304	misfdep(:,:) = jk
1305	risfdep(:,:) = zmax
1306	END WHERE
1307	END DO
1308
1309
1310	! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
1311	DO jj = 1, jpj
1312	DO ji = 1, jpi
1313	! find the minimum change option:
1314	! test bathy
1315	IF (risfdep(ji,jj) .GT. 1) THEN
1316	IF ( .NOT. ln_iscpl ) THEN
1317	zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
1318	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1319	zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &
1320	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1321
1322	IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1323	IF (zbathydiff .LE. zrisfdepdiff) THEN
1324	bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
1325	mbathy(ji,jj)= mbathy(ji,jj) + 1
1326	ELSE
1327	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1328	misfdep(ji,jj) = misfdep(ji,jj) - 1
1329	END IF
1330	ENDIF
1331	ELSE
1332	IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
1333	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
1334	misfdep(ji,jj) = misfdep(ji,jj) - 1
1335	END IF
1336	END IF
1337	END IF
1338	END DO
1339	END DO
1340
1341	! At least 2 levels for water thickness at T, U, and V point.
1342	DO jj = 1, jpj
1343	DO ji = 1, jpi
1344	! find the minimum change option:
1345	! test bathy
1346	IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1347	IF ( .NOT. ln_iscpl ) THEN
1348	zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
1349	& + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
1350	zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) &
1351	& - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
1352	IF (zbathydiff .LE. zrisfdepdiff) THEN
1353	mbathy(ji,jj) = mbathy(ji,jj) + 1
1354	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1355	ELSE
1356	misfdep(ji,jj)= misfdep(ji,jj) - 1
1357	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1358	END IF
1359	ELSE
1360	misfdep(ji,jj)= misfdep(ji,jj) - 1
1361	risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
1362	END IF
1363	ENDIF
1364	END DO
1365	END DO
1366
1367	! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)
1368	DO jj = 1, jpjm1
1369	DO ji = 1, jpim1
1370	IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1371	IF ( .NOT. ln_iscpl ) THEN
1372	zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) &
1373	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))
1374	zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
1375	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
1376	IF (zbathydiff .LE. zrisfdepdiff) THEN
1377	mbathy(ji,jj) = mbathy(ji,jj) + 1
1378	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )
1379	ELSE
1380	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
1381	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1382	END IF
1383	ELSE
1384	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
1385	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
1386	END IF
1387	ENDIF
1388	END DO
1389	END DO
1390
1391	IF( lk_mpp ) THEN
1392	zbathy(:,:) = FLOAT( misfdep(:,:) )
1393	CALL lbc_lnk( zbathy, 'T', 1. )
1394	misfdep(:,:) = INT( zbathy(:,:) )
1395	CALL lbc_lnk( risfdep, 'T', 1. )
1396	CALL lbc_lnk( bathy, 'T', 1. )
1397	zbathy(:,:) = FLOAT( mbathy(:,:) )
1398	CALL lbc_lnk( zbathy, 'T', 1. )
1399	mbathy(:,:) = INT( zbathy(:,:) )
1400	ENDIF
1401	! point V misdep(ji,jj) EQ mbathy(ji,jj+1)
1402	DO jj = 1, jpjm1
1403	DO ji = 1, jpim1
1404	IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
1405	IF ( .NOT. ln_iscpl ) THEN
1406	zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
1407	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
1408	zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) &
1409	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))
1410	IF (zbathydiff .LE. zrisfdepdiff) THEN
1411	mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
1412	bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
1413	ELSE
1414	misfdep(ji,jj) = misfdep(ji,jj) - 1
1415	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
1416	END IF
1417	ELSE
1418	misfdep(ji,jj) = misfdep(ji,jj) - 1
1419	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
1420	END IF
1421	ENDIF
1422	END DO
1423	END DO
1424
1425
1426	IF( lk_mpp ) THEN
1427	zbathy(:,:) = FLOAT( misfdep(:,:) )
1428	CALL lbc_lnk( zbathy, 'T', 1. )
1429	misfdep(:,:) = INT( zbathy(:,:) )
1430	CALL lbc_lnk( risfdep, 'T', 1. )
1431	CALL lbc_lnk( bathy, 'T', 1. )
1432	zbathy(:,:) = FLOAT( mbathy(:,:) )
1433	CALL lbc_lnk( zbathy, 'T', 1. )
1434	mbathy(:,:) = INT( zbathy(:,:) )
1435	ENDIF
1436
1437	! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)
1438	DO jj = 1, jpjm1
1439	DO ji = 1, jpim1
1440	IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1441	IF ( .NOT. ln_iscpl ) THEN
1442	zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
1443	& + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))
1444	zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
1445	& - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
1446	IF (zbathydiff .LE. zrisfdepdiff) THEN
1447	mbathy(ji,jj) = mbathy(ji,jj) + 1
1448	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
1449	ELSE
1450	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1451	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1452	END IF
1453	ELSE
1454	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
1455	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
1456	ENDIF
1457	ENDIF
1458	ENDDO
1459	ENDDO
1460
1461	IF( lk_mpp ) THEN
1462	zbathy(:,:) = FLOAT( misfdep(:,:) )
1463	CALL lbc_lnk( zbathy, 'T', 1. )
1464	misfdep(:,:) = INT( zbathy(:,:) )
1465	CALL lbc_lnk( risfdep, 'T', 1. )
1466	CALL lbc_lnk( bathy, 'T', 1. )
1467	zbathy(:,:) = FLOAT( mbathy(:,:) )
1468	CALL lbc_lnk( zbathy, 'T', 1. )
1469	mbathy(:,:) = INT( zbathy(:,:) )
1470	ENDIF
1471
1472	! point U misfdep(ji,jj) EQ bathy(ji,jj+1)
1473	DO jj = 1, jpjm1
1474	DO ji = 1, jpim1
1475	IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
1476	IF ( .NOT. ln_iscpl ) THEN
1477	zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
1478	& + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
1479	zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) &
1480	& - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))
1481	IF (zbathydiff .LE. zrisfdepdiff) THEN
1482	mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1
1483	bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
1484	ELSE
1485	misfdep(ji,jj) = misfdep(ji ,jj) - 1
1486	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
1487	END IF
1488	ELSE
1489	misfdep(ji,jj) = misfdep(ji ,jj) - 1
1490	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
1491	ENDIF
1492	ENDIF
1493	ENDDO
1494	ENDDO
1495
1496	IF( lk_mpp ) THEN
1497	zbathy(:,:) = FLOAT( misfdep(:,:) )
1498	CALL lbc_lnk( zbathy, 'T', 1. )
1499	misfdep(:,:) = INT( zbathy(:,:) )
1500	CALL lbc_lnk( risfdep, 'T', 1. )
1501	CALL lbc_lnk( bathy, 'T', 1. )
1502	zbathy(:,:) = FLOAT( mbathy(:,:) )
1503	CALL lbc_lnk( zbathy, 'T', 1. )
1504	mbathy(:,:) = INT( zbathy(:,:) )
1505	ENDIF
1506	END DO
1507	! end dig bathy/ice shelf to be compatible
1508	! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
1509	DO jl = 1,20
1510
1511	! remove single point "bay" on isf coast line in the ice shelf draft'
1512	DO jk = 2, jpk
1513	WHERE (misfdep==0) misfdep=jpk
1514	zmask=0
1515	WHERE (misfdep .LE. jk) zmask=1
1516	DO jj = 2, jpjm1
1517	DO ji = 2, jpim1
1518	IF (misfdep(ji,jj) .EQ. jk) THEN
1519	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1520	IF (ibtest .LE. 1) THEN
1521	risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
1522	IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
1523	END IF
1524	END IF
1525	END DO
1526	END DO
1527	END DO
1528	WHERE (misfdep==jpk)
1529	misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1530	END WHERE
1531	IF( lk_mpp ) THEN
1532	zbathy(:,:) = FLOAT( misfdep(:,:) )
1533	CALL lbc_lnk( zbathy, 'T', 1. )
1534	misfdep(:,:) = INT( zbathy(:,:) )
1535	CALL lbc_lnk( risfdep, 'T', 1. )
1536	CALL lbc_lnk( bathy, 'T', 1. )
1537	zbathy(:,:) = FLOAT( mbathy(:,:) )
1538	CALL lbc_lnk( zbathy, 'T', 1. )
1539	mbathy(:,:) = INT( zbathy(:,:) )
1540	ENDIF
1541
1542	! remove single point "bay" on bathy coast line beneath an ice shelf'
1543	DO jk = jpk,1,-1
1544	zmask=0
1545	WHERE (mbathy .GE. jk ) zmask=1
1546	DO jj = 2, jpjm1
1547	DO ji = 2, jpim1
1548	IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
1549	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
1550	IF (ibtest .LE. 1) THEN
1551	bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
1552	IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
1553	END IF
1554	END IF
1555	END DO
1556	END DO
1557	END DO
1558	WHERE (mbathy==0)
1559	misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
1560	END WHERE
1561	IF( lk_mpp ) THEN
1562	zbathy(:,:) = FLOAT( misfdep(:,:) )
1563	CALL lbc_lnk( zbathy, 'T', 1. )
1564	misfdep(:,:) = INT( zbathy(:,:) )
1565	CALL lbc_lnk( risfdep, 'T', 1. )
1566	CALL lbc_lnk( bathy, 'T', 1. )
1567	zbathy(:,:) = FLOAT( mbathy(:,:) )
1568	CALL lbc_lnk( zbathy, 'T', 1. )
1569	mbathy(:,:) = INT( zbathy(:,:) )
1570	ENDIF
1571
1572	! fill hole in ice shelf
1573	zmisfdep = misfdep
1574	zrisfdep = risfdep
1575	WHERE (zmisfdep .LE. 1) zmisfdep=jpk
1576	DO jj = 2, jpjm1
1577	DO ji = 2, jpim1
1578	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
1579	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
1580	IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj ) ) ibtestim1 = jpk
1581	IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj ) ) ibtestip1 = jpk
1582	IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj-1) ) ibtestjm1 = jpk
1583	IF( zmisfdep(ji,jj) .GE. mbathy(ji ,jj+1) ) ibtestjp1 = jpk
1584	ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1585	IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
1586	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
1587	END IF
1588	IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
1589	misfdep(ji,jj) = ibtest
1590	risfdep(ji,jj) = gdepw_1d(ibtest)
1591	ENDIF
1592	ENDDO
1593	ENDDO
1594
1595	IF( lk_mpp ) THEN
1596	zbathy(:,:) = FLOAT( misfdep(:,:) )
1597	CALL lbc_lnk( zbathy, 'T', 1. )
1598	misfdep(:,:) = INT( zbathy(:,:) )
1599	CALL lbc_lnk( risfdep, 'T', 1. )
1600	CALL lbc_lnk( bathy, 'T', 1. )
1601	zbathy(:,:) = FLOAT( mbathy(:,:) )
1602	CALL lbc_lnk( zbathy, 'T', 1. )
1603	mbathy(:,:) = INT( zbathy(:,:) )
1604	ENDIF
1605	!
1606	!! fill hole in bathymetry
1607	zmbathy (:,:)=mbathy (:,:)
1608	DO jj = 2, jpjm1
1609	DO ji = 2, jpim1
1610	ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )
1611	ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)
1612	IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj ) ) ibtestim1 = 0
1613	IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj ) ) ibtestip1 = 0
1614	IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj-1) ) ibtestjm1 = 0
1615	IF( zmbathy(ji,jj) .LT. misfdep(ji ,jj+1) ) ibtestjp1 = 0
1616	ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
1617	IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
1618	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
1619	END IF
1620	IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
1621	mbathy(ji,jj) = ibtest
1622	bathy(ji,jj) = gdepw_1d(ibtest+1)
1623	ENDIF
1624	END DO
1625	END DO
1626	IF( lk_mpp ) THEN
1627	zbathy(:,:) = FLOAT( misfdep(:,:) )
1628	CALL lbc_lnk( zbathy, 'T', 1. )
1629	misfdep(:,:) = INT( zbathy(:,:) )
1630	CALL lbc_lnk( risfdep, 'T', 1. )
1631	CALL lbc_lnk( bathy, 'T', 1. )
1632	zbathy(:,:) = FLOAT( mbathy(:,:) )
1633	CALL lbc_lnk( zbathy, 'T', 1. )
1634	mbathy(:,:) = INT( zbathy(:,:) )
1635	ENDIF
1636	! if not compatible after all check (ie U point water column less than 2 cells), mask U
1637	DO jj = 1, jpjm1
1638	DO ji = 1, jpim1
1639	IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1640	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
1641	END IF
1642	END DO
1643	END DO
1644	IF( lk_mpp ) THEN
1645	zbathy(:,:) = FLOAT( misfdep(:,:) )
1646	CALL lbc_lnk( zbathy, 'T', 1. )
1647	misfdep(:,:) = INT( zbathy(:,:) )
1648	CALL lbc_lnk( risfdep, 'T', 1. )
1649	CALL lbc_lnk( bathy, 'T', 1. )
1650	zbathy(:,:) = FLOAT( mbathy(:,:) )
1651	CALL lbc_lnk( zbathy, 'T', 1. )
1652	mbathy(:,:) = INT( zbathy(:,:) )
1653	ENDIF
1654	! if not compatible after all check (ie U point water column less than 2 cells), mask U
1655	DO jj = 1, jpjm1
1656	DO ji = 1, jpim1
1657	IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
1658	mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;
1659	END IF
1660	END DO
1661	END DO
1662	IF( lk_mpp ) THEN
1663	zbathy(:,:) = FLOAT( misfdep(:,:) )
1664	CALL lbc_lnk( zbathy, 'T', 1. )
1665	misfdep(:,:) = INT( zbathy(:,:) )
1666	CALL lbc_lnk( risfdep, 'T', 1. )
1667	CALL lbc_lnk( bathy, 'T', 1. )
1668	zbathy(:,:) = FLOAT( mbathy(:,:) )
1669	CALL lbc_lnk( zbathy, 'T', 1. )
1670	mbathy(:,:) = INT( zbathy(:,:) )
1671	ENDIF
1672	! if not compatible after all check (ie V point water column less than 2 cells), mask V
1673	DO jj = 1, jpjm1
1674	DO ji = 1, jpi
1675	IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1676	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
1677	END IF
1678	END DO
1679	END DO
1680	IF( lk_mpp ) THEN
1681	zbathy(:,:) = FLOAT( misfdep(:,:) )
1682	CALL lbc_lnk( zbathy, 'T', 1. )
1683	misfdep(:,:) = INT( zbathy(:,:) )
1684	CALL lbc_lnk( risfdep, 'T', 1. )
1685	CALL lbc_lnk( bathy, 'T', 1. )
1686	zbathy(:,:) = FLOAT( mbathy(:,:) )
1687	CALL lbc_lnk( zbathy, 'T', 1. )
1688	mbathy(:,:) = INT( zbathy(:,:) )
1689	ENDIF
1690	! if not compatible after all check (ie V point water column less than 2 cells), mask V
1691	DO jj = 1, jpjm1
1692	DO ji = 1, jpi
1693	IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
1694	mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
1695	END IF
1696	END DO
1697	END DO
1698	IF( lk_mpp ) THEN
1699	zbathy(:,:) = FLOAT( misfdep(:,:) )
1700	CALL lbc_lnk( zbathy, 'T', 1. )
1701	misfdep(:,:) = INT( zbathy(:,:) )
1702	CALL lbc_lnk( risfdep, 'T', 1. )
1703	CALL lbc_lnk( bathy, 'T', 1. )
1704	zbathy(:,:) = FLOAT( mbathy(:,:) )
1705	CALL lbc_lnk( zbathy, 'T', 1. )
1706	mbathy(:,:) = INT( zbathy(:,:) )
1707	ENDIF
1708	! if not compatible after all check, mask T
1709	DO jj = 1, jpj
1710	DO ji = 1, jpi
1711	IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
1712	misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;
1713	END IF
1714	END DO
1715	END DO
1716
1717	WHERE (mbathy(:,:) == 1)
1718	mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
1719	END WHERE
1720	END DO
1721	! end check compatibility ice shelf/bathy
1722	! remove very shallow ice shelf (less than ~ 10m if 75L)
1723	IF( icompt == 0 ) THEN
1724	IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'
1725	ELSE
1726	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'
1727	ENDIF
1728
1729	! compute scale factor and depth at T- and W- points
1730	DO jj = 1, jpj
1731	DO ji = 1, jpi
1732	ik = mbathy(ji,jj)
1733	IF( ik > 0 ) THEN ! ocean point only
1734	! max ocean level case
1735	IF( ik == jpkm1 ) THEN
1736	zdepwp = bathy(ji,jj)
1737	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1738	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1739	e3t_0(ji,jj,ik ) = ze3tp
1740	e3t_0(ji,jj,ik+1) = ze3tp
1741	e3w_0(ji,jj,ik ) = ze3wp
1742	e3w_0(ji,jj,ik+1) = ze3tp
1743	gdepw_0(ji,jj,ik+1) = zdepwp
1744	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1745	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1746	!
1747	ELSE ! standard case
1748	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1749	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1750	ENDIF
1751	! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1752	!gm Bug? check the gdepw_1d
1753	! ... on ik
1754	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1755	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1756	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1757	e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik )
1758	e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1)
1759	! ... on ik+1
1760	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1761	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1762	ENDIF
1763	ENDIF
1764	END DO
1765	END DO
1766	!
1767	it = 0
1768	DO jj = 1, jpj
1769	DO ji = 1, jpi
1770	ik = mbathy(ji,jj)
1771	IF( ik > 0 ) THEN ! ocean point only
1772	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1773	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1774	! test
1775	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1776	IF( zdiff <= 0._wp .AND. lwp ) THEN
1777	it = it + 1
1778	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1779	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1780	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1781	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1782	ENDIF
1783	ENDIF
1784	END DO
1785	END DO
1786	!
1787	! (ISF) Definition of e3t, u, v, w for ISF case
1788	DO jj = 1, jpj
1789	DO ji = 1, jpi
1790	ik = misfdep(ji,jj)
1791	IF( ik > 1 ) THEN ! ice shelf point only
1792	IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
1793	gdepw_0(ji,jj,ik) = risfdep(ji,jj)
1794	!gm Bug? check the gdepw_0
1795	! ... on ik
1796	gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
1797	& * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
1798	& / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
1799	e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
1800	e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1801
1802	IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
1803	e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
1804	ENDIF
1805	! ... on ik / ik-1
1806	e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik)
1807	e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1808	! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1809	gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1810	ENDIF
1811	END DO
1812	END DO
1813
1814	it = 0
1815	DO jj = 1, jpj
1816	DO ji = 1, jpi
1817	ik = misfdep(ji,jj)
1818	IF( ik > 1 ) THEN ! ice shelf point only
1819	e3tp (ji,jj) = e3t_0(ji,jj,ik )
1820	e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
1821	! test
1822	zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
1823	IF( zdiff <= 0. .AND. lwp ) THEN
1824	it = it + 1
1825	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1826	WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
1827	WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1828	WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
1829	ENDIF
1830	ENDIF
1831	END DO
1832	END DO
1833
1834	CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1835	CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1836
1837	IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')
1838
1839	END SUBROUTINE
1840
1841	SUBROUTINE zgr_sco
1842	!!----------------------------------------------------------------------
1843	!! * ROUTINE zgr_sco *
1844	!!
1845	!! ** Purpose : define the s-coordinate system
1846	!!
1847	!! ** Method : s-coordinate
1848	!! The depth of model levels is defined as the product of an
1849	!! analytical function by the local bathymetry, while the vertical
1850	!! scale factors are defined as the product of the first derivative
1851	!! of the analytical function by the bathymetry.
1852	!! (this solution save memory as depth and scale factors are not
1853	!! 3d fields)
1854	!! - Read bathymetry (in meters) at t-point and compute the
1855	!! bathymetry at u-, v-, and f-points.
1856	!! hbatu = mi( hbatt )
1857	!! hbatv = mj( hbatt )
1858	!! hbatf = mi( mj( hbatt ) )
1859	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1860	!! function and its derivative given as function.
1861	!! z_gsigt(k) = fssig (k )
1862	!! z_gsigw(k) = fssig (k-0.5)
1863	!! z_esigt(k) = fsdsig(k )
1864	!! z_esigw(k) = fsdsig(k-0.5)
1865	!! Three options for stretching are give, and they can be modified
1866	!! following the users requirements. Nevertheless, the output as
1867	!! well as the way to compute the model levels and scale factors
1868	!! must be respected in order to insure second order accuracy
1869	!! schemes.
1870	!!
1871	!! The three methods for stretching available are:
1872	!!
1873	!! s_sh94 (Song and Haidvogel 1994)
1874	!! a sinh/tanh function that allows sigma and stretched sigma
1875	!!
1876	!! s_sf12 (Siddorn and Furner 2012?)
1877	!! allows the maintenance of fixed surface and or
1878	!! bottom cell resolutions (cf. geopotential coordinates)
1879	!! within an analytically derived stretched S-coordinate framework.
1880	!!
1881	!! s_tanh (Madec et al 1996)
1882	!! a cosh/tanh function that gives stretched coordinates
1883	!!
1884	!!----------------------------------------------------------------------
1885	!
1886	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1887	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1888	INTEGER :: ios ! Local integer output status for namelist read
1889	REAL(wp) :: zrmax, ztaper ! temporary scalars
1890	REAL(wp) :: zrfact
1891	!
1892	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1893	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1894
1895	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1896	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1897	!!----------------------------------------------------------------------
1898	!
1899	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1900	!
1901	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1902	!
1903	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1904	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1905	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1906
1907	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1908	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1909	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1910	IF(lwm) WRITE ( numond, namzgr_sco )
1911
1912	IF(lwp) THEN ! control print
1913	WRITE(numout,*)
1914	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1915	WRITE(numout,*) '~~~~~~~~~~~'
1916	WRITE(numout,*) ' Namelist namzgr_sco'
1917	WRITE(numout,*) ' stretching coeffs '
1918	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1919	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1920	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1921	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1922	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1923	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1924	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1925	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1926	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1927	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1928	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1929	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1930	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1931	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1932	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1933	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1934	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1935	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1936	ENDIF
1937
1938	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1939	hifu(:,:) = rn_sbot_min
1940	hifv(:,:) = rn_sbot_min
1941	hiff(:,:) = rn_sbot_min
1942
1943	! ! set maximum ocean depth
1944	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1945
1946	DO jj = 1, jpj
1947	DO ji = 1, jpi
1948	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1949	END DO
1950	END DO
1951	! ! =============================
1952	! ! Define the envelop bathymetry (hbatt)
1953	! ! =============================
1954	! use r-value to create hybrid coordinates
1955	zenv(:,:) = bathy(:,:)
1956	!
1957	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1958	DO jj = 1, jpj
1959	DO ji = 1, jpi
1960	IF( bathy(ji,jj) == 0._wp ) THEN
1961	iip1 = MIN( ji+1, jpi )
1962	ijp1 = MIN( jj+1, jpj )
1963	iim1 = MAX( ji-1, 1 )
1964	ijm1 = MAX( jj-1, 1 )
1965	IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + &
1966	& bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1967	zenv(ji,jj) = rn_sbot_min
1968	ENDIF
1969	ENDIF
1970	END DO
1971	END DO
1972	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1973	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1974	!
1975	! smooth the bathymetry (if required)
1976	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1977	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1978	!
1979	jl = 0
1980	zrmax = 1._wp
1981	!
1982	!
1983	! set scaling factor used in reducing vertical gradients
1984	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1985	!
1986	! initialise temporary evelope depth arrays
1987	ztmpi1(:,:) = zenv(:,:)
1988	ztmpi2(:,:) = zenv(:,:)
1989	ztmpj1(:,:) = zenv(:,:)
1990	ztmpj2(:,:) = zenv(:,:)
1991	!
1992	! initialise temporary r-value arrays
1993	zri(:,:) = 1._wp
1994	zrj(:,:) = 1._wp
1995	! ! ================ !
1996	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1997	! ! ================ !
1998	jl = jl + 1
1999	zrmax = 0._wp
2000	! we set zrmax from previous r-values (zri and zrj) first
2001	! if set after current r-value calculation (as previously)
2002	! we could exit DO WHILE prematurely before checking r-value
2003	! of current zenv
2004	DO jj = 1, nlcj
2005	DO ji = 1, nlci
2006	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2007	END DO
2008	END DO
2009	zri(:,:) = 0._wp
2010	zrj(:,:) = 0._wp
2011	DO jj = 1, nlcj
2012	DO ji = 1, nlci
2013	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
2014	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
2015	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2016	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
2017	END IF
2018	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2019	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
2020	END IF
2021	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
2022	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
2023	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
2024	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
2025	END DO
2026	END DO
2027	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
2028	!
2029	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
2030	!
2031	DO jj = 1, nlcj
2032	DO ji = 1, nlci
2033	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2034	END DO
2035	END DO
2036	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
2037	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2038	! ! ================ !
2039	END DO ! End loop !
2040	! ! ================ !
2041	DO jj = 1, jpj
2042	DO ji = 1, jpi
2043	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2044	END DO
2045	END DO
2046	!
2047	! Envelope bathymetry saved in hbatt
2048	hbatt(:,:) = zenv(:,:)
2049	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2050	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2051	DO jj = 1, jpj
2052	DO ji = 1, jpi
2053	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2054	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2055	END DO
2056	END DO
2057	ENDIF
2058	!
2059	IF(lwp) THEN ! Control print
2060	WRITE(numout,*)
2061	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
2062	WRITE(numout,*)
2063	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
2064	IF( nprint == 1 ) THEN
2065	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
2066	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
2067	ENDIF
2068	ENDIF
2069
2070	! ! ==============================
2071	! ! hbatu, hbatv, hbatf fields
2072	! ! ==============================
2073	IF(lwp) THEN
2074	WRITE(numout,*)
2075	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2076	ENDIF
2077	hbatu(:,:) = rn_sbot_min
2078	hbatv(:,:) = rn_sbot_min
2079	hbatf(:,:) = rn_sbot_min
2080	DO jj = 1, jpjm1
2081	DO ji = 1, jpim1 ! NO vector opt.
2082	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
2083	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
2084	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
2085	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2086	END DO
2087	END DO
2088	!
2089	! Apply lateral boundary condition
2090	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2091	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
2092	DO jj = 1, jpj
2093	DO ji = 1, jpi
2094	IF( hbatu(ji,jj) == 0._wp ) THEN
2095	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
2096	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
2097	ENDIF
2098	END DO
2099	END DO
2100	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
2101	DO jj = 1, jpj
2102	DO ji = 1, jpi
2103	IF( hbatv(ji,jj) == 0._wp ) THEN
2104	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
2105	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
2106	ENDIF
2107	END DO
2108	END DO
2109	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
2110	DO jj = 1, jpj
2111	DO ji = 1, jpi
2112	IF( hbatf(ji,jj) == 0._wp ) THEN
2113	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
2114	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
2115	ENDIF
2116	END DO
2117	END DO
2118
2119	!!bug: key_helsinki a verifer
2120	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2121	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2122	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2123	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2124
2125	IF( nprint == 1 .AND. lwp ) THEN
2126	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
2127	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2128	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
2129	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2130	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
2131	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2132	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
2133	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2134	ENDIF
2135	!! helsinki
2136
2137	! ! =======================
2138	! ! s-ccordinate fields (gdep., e3.)
2139	! ! =======================
2140	!
2141	! non-dimensional "sigma" for model level depth at w- and t-levels
2142
2143
2144	!========================================================================
2145	! Song and Haidvogel 1994 (ln_s_sh94=T)
2146	! Siddorn and Furner 2012 (ln_sf12=T)
2147	! or tanh function (both false)
2148	!========================================================================
2149	IF ( ln_s_sh94 ) THEN
2150	CALL s_sh94()
2151	ELSE IF ( ln_s_sf12 ) THEN
2152	CALL s_sf12()
2153	ELSE
2154	CALL s_tanh()
2155	ENDIF
2156
2157	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2158	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2159	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2160	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2161	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2162	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2163	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2164
2165	fsdepw(:,:,:) = gdepw_0 (:,:,:)
2166	fsde3w(:,:,:) = gdep3w_0(:,:,:)
2167	!
2168	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
2169	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
2170	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
2171	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
2172	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
2173	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
2174	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
2175
2176	#if defined key_agrif
2177	! Ensure meaningful vertical scale factors in ghost lines/columns
2178	IF( .NOT. Agrif_Root() ) THEN
2179	!
2180	IF((nbondi == -1).OR.(nbondi == 2)) THEN
2181	e3u_0(1,:,:) = e3u_0(2,:,:)
2182	ENDIF
2183	!
2184	IF((nbondi == 1).OR.(nbondi == 2)) THEN
2185	e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
2186	ENDIF
2187	!
2188	IF((nbondj == -1).OR.(nbondj == 2)) THEN
2189	e3v_0(:,1,:) = e3v_0(:,2,:)
2190	ENDIF
2191	!
2192	IF((nbondj == 1).OR.(nbondj == 2)) THEN
2193	e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
2194	ENDIF
2195	!
2196	ENDIF
2197	#endif
2198
2199	fsdept(:,:,:) = gdept_0 (:,:,:)
2200	fsdepw(:,:,:) = gdepw_0 (:,:,:)
2201	fsde3w(:,:,:) = gdep3w_0(:,:,:)
2202	fse3t (:,:,:) = e3t_0 (:,:,:)
2203	fse3u (:,:,:) = e3u_0 (:,:,:)
2204	fse3v (:,:,:) = e3v_0 (:,:,:)
2205	fse3f (:,:,:) = e3f_0 (:,:,:)
2206	fse3w (:,:,:) = e3w_0 (:,:,:)
2207	fse3uw(:,:,:) = e3uw_0 (:,:,:)
2208	fse3vw(:,:,:) = e3vw_0 (:,:,:)
2209	!!
2210	! HYBRID :
2211	DO jj = 1, jpj
2212	DO ji = 1, jpi
2213	DO jk = 1, jpkm1
2214	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
2215	END DO
2216	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
2217	END DO
2218	END DO
2219	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
2220	& ' MAX ', MAXVAL( mbathy(:,:) )
2221
2222	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
2223	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
2224	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
2225	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
2226	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
2227	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
2228	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
2229	& ' w ', MINVAL( e3w_0 (:,:,:) )
2230
2231	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
2232	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
2233	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
2234	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
2235	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
2236	& ' w ', MAXVAL( e3w_0 (:,:,:) )
2237	ENDIF
2238	! END DO
2239	IF(lwp) THEN ! selected vertical profiles
2240	WRITE(numout,*)
2241	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2242	WRITE(numout,*) ' ~~~~~~ --------------------'
2243	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2244	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
2245	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2246	DO jj = mj0(20), mj1(20)
2247	DO ji = mi0(20), mi1(20)
2248	WRITE(numout,*)
2249	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2250	WRITE(numout,*) ' ~~~~~~ --------------------'
2251	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2252	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2253	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2254	END DO
2255	END DO
2256	DO jj = mj0(74), mj1(74)
2257	DO ji = mi0(100), mi1(100)
2258	WRITE(numout,*)
2259	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2260	WRITE(numout,*) ' ~~~~~~ --------------------'
2261	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2262	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2263	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2264	END DO
2265	END DO
2266	ENDIF
2267
2268	!================================================================================
2269	! check the coordinate makes sense
2270	!================================================================================
2271	DO ji = 1, jpi
2272	DO jj = 1, jpj
2273
2274	IF( hbatt(ji,jj) > 0._wp) THEN
2275	DO jk = 1, mbathy(ji,jj)
2276	! check coordinate is monotonically increasing
2277	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2278	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2279	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2280	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2281	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2282	CALL ctl_stop( ctmp1 )
2283	ENDIF
2284	! and check it has never gone negative
2285	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2286	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2287	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2288	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2289	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2290	CALL ctl_stop( ctmp1 )
2291	ENDIF
2292	! and check it never exceeds the total depth
2293	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2294	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2295	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2296	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2297	CALL ctl_stop( ctmp1 )
2298	ENDIF
2299	END DO
2300
2301	DO jk = 1, mbathy(ji,jj)-1
2302	! and check it never exceeds the total depth
2303	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2304	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2305	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2306	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2307	CALL ctl_stop( ctmp1 )
2308	ENDIF
2309	END DO
2310
2311	ENDIF
2312
2313	END DO
2314	END DO
2315	!
2316	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2317	!
2318	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
2319	!
2320	END SUBROUTINE zgr_sco
2321
2322	!!======================================================================
2323	SUBROUTINE s_sh94()
2324
2325	!!----------------------------------------------------------------------
2326	!! * ROUTINE s_sh94 *
2327	!!
2328	!! ** Purpose : stretch the s-coordinate system
2329	!!
2330	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
2331	!! mixed S/sigma coordinate
2332	!!
2333	!! Reference : Song and Haidvogel 1994.
2334	!!----------------------------------------------------------------------
2335	!
2336	INTEGER :: ji, jj, jk ! dummy loop argument
2337	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2338	!
2339	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2340	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2341
2342	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2343	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2344
2345	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2346	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2347	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2348	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2349
2350	DO ji = 1, jpi
2351	DO jj = 1, jpj
2352
2353	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2354	DO jk = 1, jpk
2355	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2356	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2357	END DO
2358	ELSE ! shallow water, uniform sigma
2359	DO jk = 1, jpk
2360	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
2361	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2362	END DO
2363	ENDIF
2364	!
2365	DO jk = 1, jpkm1
2366	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2367	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2368	END DO
2369	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
2370	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2371	!
2372	! Coefficients for vertical depth as the sum of e3w scale factors
2373	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2374	DO jk = 2, jpk
2375	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2376	END DO
2377	!
2378	DO jk = 1, jpk
2379	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2380	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2381	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
2382	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
2383	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
2384	END DO
2385	!
2386	END DO ! for all jj's
2387	END DO ! for all ji's
2388
2389	DO ji = 1, jpim1
2390	DO jj = 1, jpjm1
2391	DO jk = 1, jpk
2392	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
2393	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2394	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
2395	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2396	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
2397	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
2398	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2399	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
2400	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2401	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
2402	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2403	!
2404	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2405	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2406	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2407	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2408	!
2409	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2410	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2411	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2412	END DO
2413	END DO
2414	END DO
2415
2416	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2417	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2418
2419	END SUBROUTINE s_sh94
2420
2421	SUBROUTINE s_sf12
2422
2423	!!----------------------------------------------------------------------
2424	!! * ROUTINE s_sf12 *
2425	!!
2426	!! ** Purpose : stretch the s-coordinate system
2427	!!
2428	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
2429	!! mixed S/sigma/Z coordinate
2430	!!
2431	!! This method allows the maintenance of fixed surface and or
2432	!! bottom cell resolutions (cf. geopotential coordinates)
2433	!! within an analytically derived stretched S-coordinate framework.
2434	!!
2435	!!
2436	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2437	!!----------------------------------------------------------------------
2438	!
2439	INTEGER :: ji, jj, jk ! dummy loop argument
2440	REAL(wp) :: zsmth ! smoothing around critical depth
2441	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
2442	!
2443	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2444	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2445
2446	!
2447	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2448	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2449
2450	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2451	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2452	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2453	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2454
2455	DO ji = 1, jpi
2456	DO jj = 1, jpj
2457
2458	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2459
2460	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
2461	! could be changed by users but care must be taken to do so carefully
2462	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2463
2464	zzs = rn_zs / hbatt(ji,jj)
2465
2466	IF (rn_efold /= 0.0_wp) THEN
2467	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2468	ELSE
2469	zsmth = 1.0_wp
2470	ENDIF
2471
2472	DO jk = 1, jpk
2473	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2474	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2475	ENDDO
2476	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
2477	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
2478
2479	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2480
2481	DO jk = 1, jpk
2482	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2483	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2484	END DO
2485
2486	ELSE ! shallow water, z coordinates
2487
2488	DO jk = 1, jpk
2489	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2490	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2491	END DO
2492
2493	ENDIF
2494
2495	DO jk = 1, jpkm1
2496	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2497	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2498	END DO
2499	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
2500	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2501
2502	! Coefficients for vertical depth as the sum of e3w scale factors
2503	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2504	DO jk = 2, jpk
2505	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2506	END DO
2507
2508	DO jk = 1, jpk
2509	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2510	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2511	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2512	END DO
2513
2514	ENDDO ! for all jj's
2515	ENDDO ! for all ji's
2516
2517	DO ji=1,jpi-1
2518	DO jj=1,jpj-1
2519
2520	DO jk = 1, jpk
2521	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
2522	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2523	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
2524	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2525	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
2526	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
2527	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2528	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
2529	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2530	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
2531	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2532
2533	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2534	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2535	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2536	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2537	!
2538	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2539	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2540	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2541	END DO
2542
2543	ENDDO
2544	ENDDO
2545	!
2546	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2547	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2548	CALL lbc_lnk(e3w_0 ,'T',1.)
2549	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2550	!
2551	! ! =============
2552
2553	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2554	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2555
2556	END SUBROUTINE s_sf12
2557
2558	SUBROUTINE s_tanh()
2559
2560	!!----------------------------------------------------------------------
2561	!! * ROUTINE s_tanh*
2562	!!
2563	!! ** Purpose : stretch the s-coordinate system
2564	!!
2565	!! ** Method : s-coordinate stretch
2566	!!
2567	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2568	!!----------------------------------------------------------------------
2569
2570	INTEGER :: ji, jj, jk ! dummy loop argument
2571	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2572
2573	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2574	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2575
2576	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2577	CALL wrk_alloc( jpk, z_esigt, z_esigw )
2578
2579	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
2580	z_esigt = 0._wp ; z_esigw = 0._wp
2581
2582	DO jk = 1, jpk
2583	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2584	z_gsigt(jk) = -fssig( REAL(jk,wp) )
2585	END DO
2586	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
2587	!
2588	! Coefficients for vertical scale factors at w-, t- levels
2589	!!gm bug : define it from analytical function, not like juste bellow....
2590	!!gm or betteroffer the 2 possibilities....
2591	DO jk = 1, jpkm1
2592	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
2593	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2594	END DO
2595	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
2596	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2597	!
2598	! Coefficients for vertical depth as the sum of e3w scale factors
2599	z_gsi3w(1) = 0.5_wp * z_esigw(1)
2600	DO jk = 2, jpk
2601	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2602	END DO
2603	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2604	DO jk = 1, jpk
2605	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2606	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2607	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
2608	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
2609	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
2610	END DO
2611	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2612	DO jj = 1, jpj
2613	DO ji = 1, jpi
2614	DO jk = 1, jpk
2615	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2616	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2617	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2618	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2619	!
2620	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2621	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2622	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2623	END DO
2624	END DO
2625	END DO
2626
2627	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2628	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2629
2630	END SUBROUTINE s_tanh
2631
2632	FUNCTION fssig( pk ) RESULT( pf )
2633	!!----------------------------------------------------------------------
2634	!! * ROUTINE fssig *
2635	!!
2636	!! ** Purpose : provide the analytical function in s-coordinate
2637	!!
2638	!! ** Method : the function provide the non-dimensional position of
2639	!! T and W (i.e. between 0 and 1)
2640	!! T-points at integer values (between 1 and jpk)
2641	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2642	!!----------------------------------------------------------------------
2643	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2644	REAL(wp) :: pf ! sigma value
2645	!!----------------------------------------------------------------------
2646	!
2647	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2648	& - TANH( rn_thetb * rn_theta ) ) &
2649	& * ( COSH( rn_theta ) &
2650	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2651	& / ( 2._wp * SINH( rn_theta ) )
2652	!
2653	END FUNCTION fssig
2654
2655
2656	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2657	!!----------------------------------------------------------------------
2658	!! * ROUTINE fssig1 *
2659	!!
2660	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2661	!!
2662	!! ** Method : the function provides the non-dimensional position of
2663	!! T and W (i.e. between 0 and 1)
2664	!! T-points at integer values (between 1 and jpk)
2665	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2666	!!----------------------------------------------------------------------
2667	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2668	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2669	REAL(wp) :: pf1 ! sigma value
2670	!!----------------------------------------------------------------------
2671	!
2672	IF ( rn_theta == 0 ) then ! uniform sigma
2673	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2674	ELSE ! stretched sigma
2675	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2676	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2677	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2678	ENDIF
2679	!
2680	END FUNCTION fssig1
2681
2682
2683	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2684	!!----------------------------------------------------------------------
2685	!! * ROUTINE fgamma *
2686	!!
2687	!! ** Purpose : provide analytical function for the s-coordinate
2688	!!
2689	!! ** Method : the function provides the non-dimensional position of
2690	!! T and W (i.e. between 0 and 1)
2691	!! T-points at integer values (between 1 and jpk)
2692	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2693	!!
2694	!! This method allows the maintenance of fixed surface and or
2695	!! bottom cell resolutions (cf. geopotential coordinates)
2696	!! within an analytically derived stretched S-coordinate framework.
2697	!!
2698	!! Reference : Siddorn and Furner, in prep
2699	!!----------------------------------------------------------------------
2700	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2701	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2702	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2703	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2704	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2705	REAL(wp) :: za1,za2,za3 ! local variables
2706	REAL(wp) :: zn1,zn2 ! local variables
2707	REAL(wp) :: za,zb,zx ! local variables
2708	integer :: jk
2709	!!----------------------------------------------------------------------
2710	!
2711
2712	zn1 = 1./(jpk-1.)
2713	zn2 = 1. - zn1
2714
2715	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2716	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2717	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2718
2719	za = pzb - za3*(pzs-za1)-za2
2720	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2721	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2722	zx = 1.0_wp-za/2.0_wp-zb
2723
2724	DO jk = 1, jpk
2725	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2726	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2727	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2728	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2729	ENDDO
2730
2731	!
2732	END FUNCTION fgamma
2733
2734	!!======================================================================
2735	END MODULE domzgr

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

Download in other formats: