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

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

domzgr.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7509

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

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

Download in other formats: