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

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

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

Context Navigation

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5111

Last change on this file since 5111 was 5111, checked in by mathiot, 9 years ago
add some missing if ln_isfcav, test of option compatibility with ln_isfcav + small documentation
Property svn:keywords set to `Id`
File size: 130.2 KB

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

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

Download in other formats: