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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7037

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

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

Download in other formats: