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

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

domzgr.F90 in utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src – NEMO

Context Navigation

source: utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src/domzgr.F90 @ 10727

Last change on this file since 10727 was 10727, checked in by rblod, 6 years ago
new nesting tools (attempt) and brutal cleaning of DOMAINcfg, see ticket #2129
File size: 138.4 KB

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

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

Download in other formats: