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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90 @ 11658

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

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

Download in other formats: