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/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9732

Last change on this file since 9732 was 9732, checked in by mathiot, 6 years ago
replace the specific treatment of shallow ice shelves by a more general case
File size: 125.5 KB

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

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

Download in other formats: