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

Context Navigation

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7412

Last change on this file since 7412 was 7412, checked in by lovato, 8 years ago
Merge dev_NOC_CMCC_merge_2016 into branch
Property svn:keywords set to `Id`
File size: 147.5 KB

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

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

Download in other formats: