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

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

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

Context Navigation

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5862

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

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

Download in other formats: