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/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5619

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

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

Download in other formats: