Context Navigation

domzgr.F90 @ 9383

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

Download in other formats: