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 branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4724

Last change on this file since 4724 was 4724, checked in by mathiot, 10 years ago
ISF branch: add comments, fix mpp and restar issues, add test to stop if incompatible options and fix mask issue in sbcice and sbcblk.
Property svn:keywords set to `Id`
File size: 116.9 KB

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

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

Download in other formats: