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 NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src – NEMO

Context Navigation

source: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/domzgr.f90 @ 15126

Last change on this file since 15126 was 15121, checked in by dbruciaferri, 3 years ago
adding code for MEs-coordinate system v1.0 (same of PhD)
File size: 137.2 KB

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

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

Download in other formats: