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

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

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

Context Navigation

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7825

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

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

Download in other formats: