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

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

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

Context Navigation

source: branches/2015/dev_agrif_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6133

Last change on this file since 6133 was 6133, checked in by timgraham, 9 years ago
Upgraded to head of v3.6 stable XIOS2 branch
Property svn:keywords set to `Id`
File size: 130.4 KB

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

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

Download in other formats: