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

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6748

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago
GYRE hybrid parallelization
Property svn:keywords set to `Id`
File size: 141.5 KB

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

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

Download in other formats: