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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6492

Last change on this file since 6492 was 6492, checked in by mathiot, 8 years ago
Initialisation of misfdep/risfdep in case ln_zco activated + check compatibility option ln_isfcav + ln_zco/ln_sco
Property svn:keywords set to `Id`
File size: 140.0 KB

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

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

Download in other formats: