Context Navigation

domzgr.F90 @ 5428

Last change on this file since 5428 was 5428, checked in by deazer, 9 years ago
Added incode changes to allow compoarisoin with CO5 after svn keywords removed. Extracts, merges builds and runs as expected from working copy.
File size: 131.0 KB

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

Download in other formats: