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

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

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

Context Navigation

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

Last change on this file since 5200 was 5200, checked in by mathiot, 9 years ago
ISF cleaning branch: umask_i is not an interior mask, bug in definition of scale factor for bottom cell if ice shelf, remove definition of unused variables (dynhpg, ldfslp, domzgr, trasbc)
Property svn:keywords set to `Id`
File size: 133.2 KB

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

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

Download in other formats: