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

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

domzgr.F90 in branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5944

Last change on this file since 5944 was 5944, checked in by mathiot, 8 years ago
ISF: change related to reviewers comments
Property svn:keywords set to `Id`
File size: 133.0 KB

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

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

Download in other formats: