Context Navigation

domzgr.F90 @ 6934

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

Download in other formats: