Context Navigation

domzgr.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 130.9 KB

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

Download in other formats: