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

Context Navigation

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6060

Last change on this file since 6060 was 6060, checked in by timgraham, 8 years ago
Merged dev_r5836_noc2_VVL_BY_DEFAULT into branch
Property svn:keywords set to `Id`
File size: 128.0 KB

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

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

Download in other formats: