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

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

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

Download in other formats: