Context Navigation

domzgr.F90 @ 6140

Last change on this file since 6140 was 6140, checked in by timgraham, 8 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

Property svn:keywords set to Id

File size: 133.6 KB

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

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

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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6140

Download in other formats: