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

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

domzgr.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5118

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

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

Download in other formats: