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 utils/tools_dev_r11751_ENHANCE-05_SimonM-Harmonic_Analysis/DOMAINcfg/src – NEMO

Context Navigation

source: utils/tools_dev_r11751_ENHANCE-05_SimonM-Harmonic_Analysis/DOMAINcfg/src/domzgr.f90 @ 12102

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

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

Download in other formats: