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

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

domzgr.F90 in branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6404

Last change on this file since 6404 was 6404, checked in by timgraham, 8 years ago
First attempt at upgrading branch to the head of the trunk. This should include all of the simplification branch from the merge in Dec 2015.
Property svn:keywords set to `Id`
File size: 139.8 KB

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

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

Download in other formats: