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

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

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

Context Navigation

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5842

Last change on this file since 5842 was 5842, checked in by hliu, 8 years ago
Wetting and Drying update based on r:5803
Property svn:keywords set to `Id`
File size: 137.1 KB

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

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

Download in other formats: