Context Navigation

domzgr.f90 @ 15733

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

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

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

Context Navigation

source: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90 @ 15733

Download in other formats: