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

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

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

Context Navigation

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

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

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

Download in other formats: