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

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

domzgr.F90 in branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9366

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

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

Download in other formats: