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

Context Navigation

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5945

Last change on this file since 5945 was 5945, checked in by mathiot, 8 years ago
ice sheet coupling: changes based on reviewer comments
Property svn:keywords set to `Id`
File size: 134.9 KB

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

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

Download in other formats: