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

Context Navigation

source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6986

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

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

Download in other formats: