Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domzgr.F90 @ 9529

Last change on this file since 9529 was 9529, checked in by jchanut, 6 years ago
ztilde stability update
Property svn:keywords set to `Id`
File size: 130.7 KB

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

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

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

Context Navigation

source: branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9529

Download in other formats: