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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6152

Last change on this file since 6152 was 6152, checked in by acc, 8 years ago
Add wetting and drying option from dev_r5803_NOC_WAD branch. Logically isolated code changes in domvvl.F90, domzgr.F90, dynhpg.F90, dynspg_ts.F90, sshwzv.F90 and nemogcm.F90. New module wet_dry.F90 in DYN. Fully SETTE tested with code deactivated (ln_wad=.false.). No test case yet available to justify activating option (still under development)
Property svn:keywords set to `Id`
File size: 139.6 KB

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

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

Download in other formats: