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

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6624

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

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

Download in other formats: