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

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

domzgr.F90 in branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

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

Last change on this file since 5779 was 5779, checked in by mathiot, 9 years ago
ISF coupling branch: correct some compilation issues, remove code related to MISOMIP/ISOMIP+ and polishing
Property svn:keywords set to `Id`
File size: 132.6 KB

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

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

Download in other formats: