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 @ 9513

Last change on this file since 9513 was 9513, checked in by mathiot, 6 years ago
Add option to detect and remove subglacial lake (do not affect closed sea option)
File size: 125.2 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
7	!! ! 1997-07 (G. Madec) lbc_lnk call
8	!! ! 1997-04 (J.-O. Beismann)
9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
19	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
20	!!----------------------------------------------------------------------
21
22	!!----------------------------------------------------------------------
23	!! dom_zgr : defined the ocean vertical coordinate system
24	!! zgr_bat : bathymetry fields (levels and meters)
25	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
26	!! zgr_bat_ctl : check the bathymetry files
27	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
28	!! zgr_z : reference z-coordinate
29	!! zgr_zco : z-coordinate
30	!! zgr_zps : z-coordinate with partial steps
31	!! zgr_sco : s-coordinate
32	!! fssig : tanh stretch function
33	!! fssig1 : Song and Haidvogel 1994 stretch function
34	!! fgamma : Siddorn and Furner 2012 stretching function
35	!!---------------------------------------------------------------------
36	USE oce ! ocean variables
37	USE dom_oce ! ocean domain
38	USE closea ! closed seas
39	USE c1d ! 1D vertical configuration
40	USE in_out_manager ! I/O manager
41	USE iom ! I/O library
42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
43	USE lib_mpp ! distributed memory computing library
44	USE wrk_nemo ! Memory allocation
45	USE timing ! Timing
46
47	IMPLICIT NONE
48	PRIVATE
49
50	PUBLIC dom_zgr ! called by dom_init.F90
51
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 in shallow water
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_isfhw_deep = 1.e-3 , rn_isfhw_shallow = 1.0e-3 ! threshold to define grounding line in deep/shallow water
1301	REAL(wp) :: rn_isfshallow = 0.0_wp ! threshold to define shallow ice shelf cavity
1302	REAL(wp) :: rn_zisfmax = 6000.0_wp ! maximun meter of ice we are allowed to dig to assure connectivity
1303	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)
1304	!!---------------------------------------------------------------------
1305	NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
1306	& rn_isfdep_min, rn_isfhw_deep, rn_isfhw_shallow, rn_isfshallow, &
1307	& ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
1308	!
1309	IF( nn_timing == 1 ) CALL timing_start('zgr_isf')
1310	!
1311	CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep)
1312	CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1313	!
1314	!
1315	REWIND( numnam_ref ) ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition
1316	READ ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901)
1317	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp )
1318
1319	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition
1320	READ ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 )
1321	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp )
1322	IF(lwm) WRITE ( numond, namzgr_isf )
1323	!
1324	IF (lwp) THEN
1325	WRITE(numout,*) ' nn_kisfmax = ',nn_kisfmax !
1326	WRITE(numout,*) ' rn_zisfmax = ',rn_zisfmax !
1327	WRITE(numout,*) ' rn_isfdep_min = ',rn_isfdep_min !
1328	WRITE(numout,*) ' rn_isfhw_deep = ',rn_isfhw_deep !
1329	WRITE(numout,*) ' rn_isfhw_shallow = ',rn_isfhw_shallow !
1330	WRITE(numout,*) ' rn_isfshallow = ',rn_isfshallow !
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
1342	WHERE( bathy(:,:) - risfdep(:,:) < rn_isfhw_deep .AND. (bathy(:,:) > rn_isfshallow) )
1343	risfdep(:,:)=0.0 ; misfdep(:,:) = 1
1344	bathy (:,:)=0.0 ; mbathy (:,:) = 0
1345	END WHERE
1346	!
1347	! 1.2 ensure a minimum thickness for iceshelf cavity in shallow water
1348	WHERE( (bathy(:,:) <= rn_isfshallow) .AND. (bathy(:,:) - risfdep(:,:) < rn_isfhw_shallow) )
1349	risfdep(:,:)=0.0 ; misfdep(:,:) = 1
1350	bathy (:,:)=0.0 ; mbathy (:,:) = 0
1351	END WHERE
1352	!
1353	! 1.3 Remove channels and single point 'bay'.
1354	! channel could be created if connectivity is enforced later.
1355	IF (ln_isfchannel) THEN
1356	zmask(:,:)=1._wp
1357	WHERE (mbathy(:,:)==0.0)
1358	zmask(:,:)=0._wp
1359	END WHERE
1360	DO jj = 2, jpjm1
1361	DO ji = 2, jpim1
1362	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1363	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = connexion
1364	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = connexion
1365	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = connexion
1366	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = connexion
1367	! zonal channel and single bay
1368	IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN
1369	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1370	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1371	END IF
1372	! meridionnal channel and single bay
1373	IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN
1374	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1375	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1376	END IF
1377	END IF
1378	END DO
1379	END DO
1380	! ensure halo correct
1381	IF( lk_mpp ) THEN
1382	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1383	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1384	CALL lbc_lnk( risfdep, 'T', 1. )
1385	CALL lbc_lnk( bathy , 'T', 1. )
1386	ENDIF
1387	END IF
1388	!
1389	! 1.4 Closed single isolated point
1390	zmask(:,:)=1._wp
1391	WHERE (mbathy(:,:)==0.0)
1392	zmask(:,:)=0._wp
1393	END WHERE
1394	DO jj = 2, jpjm1
1395	DO ji = 2, jpim1
1396	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1397	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = connexion
1398	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = connexion
1399	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = connexion
1400	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = connexion
1401	! remove single point
1402	IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN
1403	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1404	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1405	END IF
1406	END IF
1407	END DO
1408	END DO
1409	! ensure halo correct
1410	IF( lk_mpp ) THEN
1411	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1412	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1413	CALL lbc_lnk( risfdep, 'T', 1. )
1414	CALL lbc_lnk( bathy , 'T', 1. )
1415	ENDIF
1416	!
1417	! 2 Compute misfdep for ocean points (i.e. first wet level)
1418	! find the first ocean level such that the first level thickness
1419	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1420	! e3t_0 is the reference level thickness
1421	!
1422	! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf
1423	WHERE( risfdep(:,:) == 0.0 ) ; misfdep(:,:) = 1 ! open water or grounded : set misfdep to 1
1424	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
1425	END WHERE
1426	!
1427	DO jk = 2, jpkm1
1428	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1429	WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
1430	END DO
1431	!
1432	! 2.1 be sure misfdep not below mbathy
1433	! warning because of condition 4 we could have 'wet cell with misfdep below mbathy
1434	! risfdep of these cells will be fix later on (see 5)
1435	WHERE( misfdep > mbathy .AND. mbathy > 0 ) misfdep=mbathy
1436	!
1437	! 3.0 Assure 2 wet cells in the column at T point and along the edge.
1438	! first do T, then loop 4 times on U-1, U, V-1, V and then at T point
1439	! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
1440	!
1441	! find the deepest isfdep level that fit the 2 wet cell on the water column
1442	! on all the sides (still need 4 pass)
1443	DO jn = 1, 4
1444	zmisfdep = misfdep
1445	DO jj = 2, jpjm1
1446	DO ji = 2, jpim1
1447	! ISF cell only
1448	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1449	zmbathyij = mbathy(ji ,jj)
1450	! set ground edge value to jpk to skip it later on
1451	zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge
1452	zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge
1453	zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge
1454	zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge
1455	! shallowest bathy level
1456	zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1)
1457	! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the
1458	! water column
1459	jk = MIN(misfdep(ji,jj),zmbathyij-1)
1460	! update isf draft if needed
1461	IF ( jk < misfdep(ji,jj) ) THEN
1462	zmisfdep(ji,jj) = jk
1463	risfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1464	ELSE
1465	! isf depth not modify, nothing to do
1466	END IF
1467	ENDIF
1468	END DO
1469	END DO
1470	misfdep=zmisfdep
1471	! ensure halo correct before new pass
1472	IF( lk_mpp ) THEN
1473	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1474	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1475	CALL lbc_lnk( risfdep, 'T', 1. )
1476	CALL lbc_lnk( bathy , 'T', 1. )
1477	ENDIF
1478	END DO ! jn
1479	!
1480	! 3.1 condition block to assure connectivity every where beneath an ice shelf
1481	IF (ln_isfconnect) THEN
1482	zmask(:,:)=1.0
1483	WHERE (mbathy==0)
1484	zmask(:,:)=jpk
1485	END WHERE
1486
1487	zmbathy = mbathy
1488	zmisfdep = misfdep
1489	zrisfdep = risfdep
1490	WHERE (mbathy == 0) zmbathy=jpk
1491	DO jj = 2, jpjm1
1492	DO ji = 2, jpim1
1493	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1494	! what it should be
1495	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = should be connected
1496	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = should be connected
1497	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = should be connected
1498	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = should be connected
1499	! what it is
1500	umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk
1501	IF (misfdep(ji,jj) > zmbathy(ji+1,jj )) umskip1_r=1.0 ! 1 = no effective connection
1502	IF (misfdep(ji,jj) > zmbathy(ji-1,jj )) umskim1_r=1.0
1503	IF (misfdep(ji,jj) > zmbathy(ji ,jj+1)) vmskjp1_r=1.0
1504	IF (misfdep(ji,jj) > zmbathy(ji ,jj-1)) vmskjm1_r=1.0
1505	! defining level need for connectivity
1506	jk=NINT(MIN(zmbathy(ji+1,jj ) * umskip1_r * umskip1, &
1507	& zmbathy(ji-1,jj ) * umskim1_r * umskim1, &
1508	& zmbathy(ji ,jj+1) * vmskjp1_r * vmskjp1, &
1509	& zmbathy(ji ,jj-1) * vmskjm1_r * vmskjm1, &
1510	& jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on
1511	! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
1512	zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
1513	!
1514	! update isf draft if needed (need to be done here because next test is in meter and level)
1515	IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN
1516	jk = zmisfdep(ji,jj)
1517	zrisfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1518	END IF
1519
1520	! sanity check
1521	! check if we dig more than nn_kisfmax level or reach the surface
1522	! check if we dig more than rn_zisfmax meter
1523	! => if this is the case, undo what has been done before
1524	IF ( (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
1525	& .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj) )) ) THEN
1526	zmisfdep(ji,jj)=misfdep(ji,jj)
1527	zrisfdep(ji,jj)=risfdep(ji,jj)
1528	END IF
1529	END IF
1530	END DO
1531	END DO
1532	misfdep=zmisfdep
1533	risfdep=zrisfdep
1534	! ensure halo correct
1535	IF( lk_mpp ) THEN
1536	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1537	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1538	CALL lbc_lnk( risfdep, 'T', 1. )
1539	CALL lbc_lnk( bathy , 'T', 1. )
1540	ENDIF
1541	END IF
1542	!
1543	IF (ln_isfcheminey) THEN
1544	! 3.2 fill hole in ice shelf (ie cell with no velocity point)
1545	! => misfdep = MIN(misfdep at U, U-1, V, V-1)
1546	! risfdep = gdepw(misfdep) (ie full cell)
1547	zmisfdep = misfdep
1548	WHERE (mbathy == 0) zmisfdep=jpk ! grounded
1549	DO jj = 2, jpjm1
1550	DO ji = 2, jpim1
1551	ibtest = zmisfdep(ji ,jj )
1552	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
1553	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
1554	! case unconected wet cell (T point) where isfdraft(ii,jj) deeper
1555	! than bathy U-1/U+1/V-1/V+1
1556	IF( zmisfdep(ji,jj) > mbathy(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1557	IF( zmisfdep(ji,jj) > mbathy(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1558	IF( zmisfdep(ji,jj) > mbathy(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1559	IF( zmisfdep(ji,jj) > mbathy(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1560	! case unconected wet cell (T point) where bathy(ii,jj) shallower
1561	! than isfdraft U-1/U+1/V-1/V+1
1562	IF( mbathy(ji,jj) < zmisfdep(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1563	IF( mbathy(ji,jj) < zmisfdep(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1564	IF( mbathy(ji,jj) < zmisfdep(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1565	IF( mbathy(ji,jj) < zmisfdep(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1566	! if surround in fully mask => mask cell (1)
1567	! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)
1568	! if misfdep is not changed, nothing is done (the MAX is doing the
1569	! job)
1570	ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
1571	IF (misfdep(ji,jj) < ibtest) THEN
1572	misfdep(ji,jj) = ibtest
1573	risfdep(ji,jj) = gdepw_1d(ibtest)
1574	END IF
1575	ENDDO
1576	ENDDO
1577	WHERE( misfdep == jpk) ! ground case (1)
1578	mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp
1579	END WHERE
1580	!
1581	! ensure halo correct
1582	IF( lk_mpp ) THEN
1583	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1584	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1585	CALL lbc_lnk( risfdep, 'T', 1. )
1586	CALL lbc_lnk( bathy , 'T', 1. )
1587	END IF
1588	END IF
1589	!
1590	CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep)
1591	CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy)
1592	!
1593	IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')
1594
1595	END SUBROUTINE
1596
1597	SUBROUTINE zgr_sco
1598	!!----------------------------------------------------------------------
1599	!! * ROUTINE zgr_sco *
1600	!!
1601	!! ** Purpose : define the s-coordinate system
1602	!!
1603	!! ** Method : s-coordinate
1604	!! The depth of model levels is defined as the product of an
1605	!! analytical function by the local bathymetry, while the vertical
1606	!! scale factors are defined as the product of the first derivative
1607	!! of the analytical function by the bathymetry.
1608	!! (this solution save memory as depth and scale factors are not
1609	!! 3d fields)
1610	!! - Read bathymetry (in meters) at t-point and compute the
1611	!! bathymetry at u-, v-, and f-points.
1612	!! hbatu = mi( hbatt )
1613	!! hbatv = mj( hbatt )
1614	!! hbatf = mi( mj( hbatt ) )
1615	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1616	!! function and its derivative given as function.
1617	!! z_gsigt(k) = fssig (k )
1618	!! z_gsigw(k) = fssig (k-0.5)
1619	!! z_esigt(k) = fsdsig(k )
1620	!! z_esigw(k) = fsdsig(k-0.5)
1621	!! Three options for stretching are give, and they can be modified
1622	!! following the users requirements. Nevertheless, the output as
1623	!! well as the way to compute the model levels and scale factors
1624	!! must be respected in order to insure second order accuracy
1625	!! schemes.
1626	!!
1627	!! The three methods for stretching available are:
1628	!!
1629	!! s_sh94 (Song and Haidvogel 1994)
1630	!! a sinh/tanh function that allows sigma and stretched sigma
1631	!!
1632	!! s_sf12 (Siddorn and Furner 2012?)
1633	!! allows the maintenance of fixed surface and or
1634	!! bottom cell resolutions (cf. geopotential coordinates)
1635	!! within an analytically derived stretched S-coordinate framework.
1636	!!
1637	!! s_tanh (Madec et al 1996)
1638	!! a cosh/tanh function that gives stretched coordinates
1639	!!
1640	!!----------------------------------------------------------------------
1641	!
1642	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1643	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1644	INTEGER :: ios ! Local integer output status for namelist read
1645	REAL(wp) :: zrmax, ztaper ! temporary scalars
1646	REAL(wp) :: zrfact
1647	!
1648	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1649	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1650
1651	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1652	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1653	!!----------------------------------------------------------------------
1654	!
1655	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1656	!
1657	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1658	!
1659	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1660	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1661	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1662
1663	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1664	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1665	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1666	IF(lwm) WRITE ( numond, namzgr_sco )
1667
1668	IF(lwp) THEN ! control print
1669	WRITE(numout,*)
1670	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1671	WRITE(numout,*) '~~~~~~~~~~~'
1672	WRITE(numout,*) ' Namelist namzgr_sco'
1673	WRITE(numout,*) ' stretching coeffs '
1674	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1675	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1676	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1677	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1678	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1679	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1680	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1681	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1682	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1683	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1684	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1685	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1686	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1687	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1688	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1689	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1690	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1691	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1692	ENDIF
1693
1694	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1695	hifu(:,:) = rn_sbot_min
1696	hifv(:,:) = rn_sbot_min
1697	hiff(:,:) = rn_sbot_min
1698
1699	! ! set maximum ocean depth
1700	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1701
1702	DO jj = 1, jpj
1703	DO ji = 1, jpi
1704	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1705	END DO
1706	END DO
1707	! ! =============================
1708	! ! Define the envelop bathymetry (hbatt)
1709	! ! =============================
1710	! use r-value to create hybrid coordinates
1711	zenv(:,:) = bathy(:,:)
1712	!
1713	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1714	DO jj = 1, jpj
1715	DO ji = 1, jpi
1716	IF( bathy(ji,jj) == 0._wp ) THEN
1717	iip1 = MIN( ji+1, jpi )
1718	ijp1 = MIN( jj+1, jpj )
1719	iim1 = MAX( ji-1, 1 )
1720	ijm1 = MAX( jj-1, 1 )
1721	IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
1722	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
1723	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) THEN
1724	zenv(ji,jj) = rn_sbot_min
1725	ENDIF
1726	ENDIF
1727	END DO
1728	END DO
1729	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1730	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1731	!
1732	! smooth the bathymetry (if required)
1733	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1734	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1735	!
1736	jl = 0
1737	zrmax = 1._wp
1738	!
1739	!
1740	! set scaling factor used in reducing vertical gradients
1741	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1742	!
1743	! initialise temporary evelope depth arrays
1744	ztmpi1(:,:) = zenv(:,:)
1745	ztmpi2(:,:) = zenv(:,:)
1746	ztmpj1(:,:) = zenv(:,:)
1747	ztmpj2(:,:) = zenv(:,:)
1748	!
1749	! initialise temporary r-value arrays
1750	zri(:,:) = 1._wp
1751	zrj(:,:) = 1._wp
1752	! ! ================ !
1753	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1754	! ! ================ !
1755	jl = jl + 1
1756	zrmax = 0._wp
1757	! we set zrmax from previous r-values (zri and zrj) first
1758	! if set after current r-value calculation (as previously)
1759	! we could exit DO WHILE prematurely before checking r-value
1760	! of current zenv
1761	DO jj = 1, nlcj
1762	DO ji = 1, nlci
1763	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1764	END DO
1765	END DO
1766	zri(:,:) = 0._wp
1767	zrj(:,:) = 0._wp
1768	DO jj = 1, nlcj
1769	DO ji = 1, nlci
1770	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1771	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1772	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1773	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1774	END IF
1775	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1776	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1777	END IF
1778	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
1779	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
1780	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
1781	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
1782	END DO
1783	END DO
1784	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1785	!
1786	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
1787	!
1788	DO jj = 1, nlcj
1789	DO ji = 1, nlci
1790	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1791	END DO
1792	END DO
1793	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1794	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1795	! ! ================ !
1796	END DO ! End loop !
1797	! ! ================ !
1798	DO jj = 1, jpj
1799	DO ji = 1, jpi
1800	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1801	END DO
1802	END DO
1803	!
1804	! Envelope bathymetry saved in hbatt
1805	hbatt(:,:) = zenv(:,:)
1806	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1807	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1808	DO jj = 1, jpj
1809	DO ji = 1, jpi
1810	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1811	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1812	END DO
1813	END DO
1814	ENDIF
1815	!
1816	IF(lwp) THEN ! Control print
1817	WRITE(numout,*)
1818	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1819	WRITE(numout,*)
1820	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1821	IF( nprint == 1 ) THEN
1822	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1823	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1824	ENDIF
1825	ENDIF
1826
1827	! ! ==============================
1828	! ! hbatu, hbatv, hbatf fields
1829	! ! ==============================
1830	IF(lwp) THEN
1831	WRITE(numout,*)
1832	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1833	ENDIF
1834	hbatu(:,:) = rn_sbot_min
1835	hbatv(:,:) = rn_sbot_min
1836	hbatf(:,:) = rn_sbot_min
1837	DO jj = 1, jpjm1
1838	DO ji = 1, jpim1 ! NO vector opt.
1839	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1840	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1841	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1842	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1843	END DO
1844	END DO
1845	!
1846	! Apply lateral boundary condition
1847	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1848	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
1849	DO jj = 1, jpj
1850	DO ji = 1, jpi
1851	IF( hbatu(ji,jj) == 0._wp ) THEN
1852	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1853	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1854	ENDIF
1855	END DO
1856	END DO
1857	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
1858	DO jj = 1, jpj
1859	DO ji = 1, jpi
1860	IF( hbatv(ji,jj) == 0._wp ) THEN
1861	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1862	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1863	ENDIF
1864	END DO
1865	END DO
1866	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
1867	DO jj = 1, jpj
1868	DO ji = 1, jpi
1869	IF( hbatf(ji,jj) == 0._wp ) THEN
1870	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1871	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1872	ENDIF
1873	END DO
1874	END DO
1875
1876	!!bug: key_helsinki a verifer
1877	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1878	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1879	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1880	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1881
1882	IF( nprint == 1 .AND. lwp ) THEN
1883	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1884	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1885	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1886	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1887	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1888	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1889	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1890	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1891	ENDIF
1892	!! helsinki
1893
1894	! ! =======================
1895	! ! s-ccordinate fields (gdep., e3.)
1896	! ! =======================
1897	!
1898	! non-dimensional "sigma" for model level depth at w- and t-levels
1899
1900
1901	!========================================================================
1902	! Song and Haidvogel 1994 (ln_s_sh94=T)
1903	! Siddorn and Furner 2012 (ln_sf12=T)
1904	! or tanh function (both false)
1905	!========================================================================
1906	IF ( ln_s_sh94 ) THEN
1907	CALL s_sh94()
1908	ELSE IF ( ln_s_sf12 ) THEN
1909	CALL s_sf12()
1910	ELSE
1911	CALL s_tanh()
1912	ENDIF
1913
1914	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
1915	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
1916	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
1917	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
1918	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
1919	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
1920	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1921
1922	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1923	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1924	!
1925	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
1926	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
1927	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
1928	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
1929	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
1930	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
1931	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
1932
1933	#if defined key_agrif
1934	! Ensure meaningful vertical scale factors in ghost lines/columns
1935	IF( .NOT. Agrif_Root() ) THEN
1936	!
1937	IF((nbondi == -1).OR.(nbondi == 2)) THEN
1938	e3u_0(1,:,:) = e3u_0(2,:,:)
1939	ENDIF
1940	!
1941	IF((nbondi == 1).OR.(nbondi == 2)) THEN
1942	e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
1943	ENDIF
1944	!
1945	IF((nbondj == -1).OR.(nbondj == 2)) THEN
1946	e3v_0(:,1,:) = e3v_0(:,2,:)
1947	ENDIF
1948	!
1949	IF((nbondj == 1).OR.(nbondj == 2)) THEN
1950	e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
1951	ENDIF
1952	!
1953	ENDIF
1954	#endif
1955
1956	fsdept(:,:,:) = gdept_0 (:,:,:)
1957	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1958	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1959	fse3t (:,:,:) = e3t_0 (:,:,:)
1960	fse3u (:,:,:) = e3u_0 (:,:,:)
1961	fse3v (:,:,:) = e3v_0 (:,:,:)
1962	fse3f (:,:,:) = e3f_0 (:,:,:)
1963	fse3w (:,:,:) = e3w_0 (:,:,:)
1964	fse3uw(:,:,:) = e3uw_0 (:,:,:)
1965	fse3vw(:,:,:) = e3vw_0 (:,:,:)
1966	!!
1967	! HYBRID :
1968	DO jj = 1, jpj
1969	DO ji = 1, jpi
1970	DO jk = 1, jpkm1
1971	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1972	END DO
1973	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
1974	END DO
1975	END DO
1976	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1977	& ' MAX ', MAXVAL( mbathy(:,:) )
1978
1979	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1980	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1981	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
1982	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
1983	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
1984	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
1985	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
1986	& ' w ', MINVAL( e3w_0 (:,:,:) )
1987
1988	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
1989	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
1990	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
1991	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
1992	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
1993	& ' w ', MAXVAL( e3w_0 (:,:,:) )
1994	ENDIF
1995	! END DO
1996	IF(lwp) THEN ! selected vertical profiles
1997	WRITE(numout,*)
1998	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1999	WRITE(numout,*) ' ~~~~~~ --------------------'
2000	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2001	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
2002	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2003	DO jj = mj0(20), mj1(20)
2004	DO ji = mi0(20), mi1(20)
2005	WRITE(numout,*)
2006	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2007	WRITE(numout,*) ' ~~~~~~ --------------------'
2008	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2009	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2010	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2011	END DO
2012	END DO
2013	DO jj = mj0(74), mj1(74)
2014	DO ji = mi0(100), mi1(100)
2015	WRITE(numout,*)
2016	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,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	ENDIF
2024
2025	!================================================================================
2026	! check the coordinate makes sense
2027	!================================================================================
2028	DO ji = 1, jpi
2029	DO jj = 1, jpj
2030
2031	IF( hbatt(ji,jj) > 0._wp) THEN
2032	DO jk = 1, mbathy(ji,jj)
2033	! check coordinate is monotonically increasing
2034	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2035	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2036	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2037	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2038	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2039	CALL ctl_stop( ctmp1 )
2040	ENDIF
2041	! and check it has never gone negative
2042	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2043	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2044	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2045	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2046	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2047	CALL ctl_stop( ctmp1 )
2048	ENDIF
2049	! and check it never exceeds the total depth
2050	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2051	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2052	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2053	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2054	CALL ctl_stop( ctmp1 )
2055	ENDIF
2056	END DO
2057
2058	DO jk = 1, mbathy(ji,jj)-1
2059	! and check it never exceeds the total depth
2060	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2061	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2062	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2063	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2064	CALL ctl_stop( ctmp1 )
2065	ENDIF
2066	END DO
2067
2068	ENDIF
2069
2070	END DO
2071	END DO
2072	!
2073	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2074	!
2075	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
2076	!
2077	END SUBROUTINE zgr_sco
2078
2079	!!======================================================================
2080	SUBROUTINE s_sh94()
2081
2082	!!----------------------------------------------------------------------
2083	!! * ROUTINE s_sh94 *
2084	!!
2085	!! ** Purpose : stretch the s-coordinate system
2086	!!
2087	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
2088	!! mixed S/sigma coordinate
2089	!!
2090	!! Reference : Song and Haidvogel 1994.
2091	!!----------------------------------------------------------------------
2092	!
2093	INTEGER :: ji, jj, jk ! dummy loop argument
2094	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2095	!
2096	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2097	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2098
2099	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2100	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2101
2102	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2103	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2104	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2105	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2106
2107	DO ji = 1, jpi
2108	DO jj = 1, jpj
2109
2110	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2111	DO jk = 1, jpk
2112	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2113	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2114	END DO
2115	ELSE ! shallow water, uniform sigma
2116	DO jk = 1, jpk
2117	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
2118	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2119	END DO
2120	ENDIF
2121	!
2122	DO jk = 1, jpkm1
2123	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2124	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2125	END DO
2126	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
2127	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2128	!
2129	! Coefficients for vertical depth as the sum of e3w scale factors
2130	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2131	DO jk = 2, jpk
2132	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2133	END DO
2134	!
2135	DO jk = 1, jpk
2136	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2137	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2138	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
2139	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
2140	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
2141	END DO
2142	!
2143	END DO ! for all jj's
2144	END DO ! for all ji's
2145
2146	DO ji = 1, jpim1
2147	DO jj = 1, jpjm1
2148	DO jk = 1, jpk
2149	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
2150	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2151	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
2152	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2153	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
2154	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
2155	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2156	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
2157	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2158	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
2159	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2160	!
2161	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2162	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2163	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2164	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2165	!
2166	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2167	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2168	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2169	END DO
2170	END DO
2171	END DO
2172
2173	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2174	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2175
2176	END SUBROUTINE s_sh94
2177
2178	SUBROUTINE s_sf12
2179
2180	!!----------------------------------------------------------------------
2181	!! * ROUTINE s_sf12 *
2182	!!
2183	!! ** Purpose : stretch the s-coordinate system
2184	!!
2185	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
2186	!! mixed S/sigma/Z coordinate
2187	!!
2188	!! This method allows the maintenance of fixed surface and or
2189	!! bottom cell resolutions (cf. geopotential coordinates)
2190	!! within an analytically derived stretched S-coordinate framework.
2191	!!
2192	!!
2193	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2194	!!----------------------------------------------------------------------
2195	!
2196	INTEGER :: ji, jj, jk ! dummy loop argument
2197	REAL(wp) :: zsmth ! smoothing around critical depth
2198	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
2199	!
2200	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2201	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2202
2203	!
2204	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2205	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2206
2207	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2208	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2209	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2210	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2211
2212	DO ji = 1, jpi
2213	DO jj = 1, jpj
2214
2215	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2216
2217	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
2218	! could be changed by users but care must be taken to do so carefully
2219	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2220
2221	zzs = rn_zs / hbatt(ji,jj)
2222
2223	IF (rn_efold /= 0.0_wp) THEN
2224	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2225	ELSE
2226	zsmth = 1.0_wp
2227	ENDIF
2228
2229	DO jk = 1, jpk
2230	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2231	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2232	ENDDO
2233	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
2234	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
2235
2236	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2237
2238	DO jk = 1, jpk
2239	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2240	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2241	END DO
2242
2243	ELSE ! shallow water, z coordinates
2244
2245	DO jk = 1, jpk
2246	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2247	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2248	END DO
2249
2250	ENDIF
2251
2252	DO jk = 1, jpkm1
2253	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2254	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2255	END DO
2256	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
2257	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2258
2259	! Coefficients for vertical depth as the sum of e3w scale factors
2260	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2261	DO jk = 2, jpk
2262	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2263	END DO
2264
2265	DO jk = 1, jpk
2266	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2267	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2268	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2269	END DO
2270
2271	ENDDO ! for all jj's
2272	ENDDO ! for all ji's
2273
2274	DO ji=1,jpi-1
2275	DO jj=1,jpj-1
2276
2277	DO jk = 1, jpk
2278	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
2279	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2280	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
2281	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2282	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
2283	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
2284	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2285	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
2286	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2287	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
2288	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2289
2290	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2291	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2292	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2293	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2294	!
2295	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2296	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2297	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2298	END DO
2299
2300	ENDDO
2301	ENDDO
2302	!
2303	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2304	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2305	CALL lbc_lnk(e3w_0 ,'T',1.)
2306	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2307	!
2308	! ! =============
2309
2310	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2311	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2312
2313	END SUBROUTINE s_sf12
2314
2315	SUBROUTINE s_tanh()
2316
2317	!!----------------------------------------------------------------------
2318	!! * ROUTINE s_tanh*
2319	!!
2320	!! ** Purpose : stretch the s-coordinate system
2321	!!
2322	!! ** Method : s-coordinate stretch
2323	!!
2324	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2325	!!----------------------------------------------------------------------
2326
2327	INTEGER :: ji, jj, jk ! dummy loop argument
2328	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2329
2330	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2331	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2332
2333	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2334	CALL wrk_alloc( jpk, z_esigt, z_esigw )
2335
2336	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
2337	z_esigt = 0._wp ; z_esigw = 0._wp
2338
2339	DO jk = 1, jpk
2340	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2341	z_gsigt(jk) = -fssig( REAL(jk,wp) )
2342	END DO
2343	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
2344	!
2345	! Coefficients for vertical scale factors at w-, t- levels
2346	!!gm bug : define it from analytical function, not like juste bellow....
2347	!!gm or betteroffer the 2 possibilities....
2348	DO jk = 1, jpkm1
2349	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
2350	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2351	END DO
2352	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
2353	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2354	!
2355	! Coefficients for vertical depth as the sum of e3w scale factors
2356	z_gsi3w(1) = 0.5_wp * z_esigw(1)
2357	DO jk = 2, jpk
2358	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2359	END DO
2360	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2361	DO jk = 1, jpk
2362	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2363	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2364	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
2365	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
2366	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
2367	END DO
2368	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2369	DO jj = 1, jpj
2370	DO ji = 1, jpi
2371	DO jk = 1, jpk
2372	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2373	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2374	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2375	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2376	!
2377	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2378	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2379	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2380	END DO
2381	END DO
2382	END DO
2383
2384	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2385	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2386
2387	END SUBROUTINE s_tanh
2388
2389	FUNCTION fssig( pk ) RESULT( pf )
2390	!!----------------------------------------------------------------------
2391	!! * ROUTINE fssig *
2392	!!
2393	!! ** Purpose : provide the analytical function in s-coordinate
2394	!!
2395	!! ** Method : the function provide the non-dimensional position of
2396	!! T and W (i.e. between 0 and 1)
2397	!! T-points at integer values (between 1 and jpk)
2398	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2399	!!----------------------------------------------------------------------
2400	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2401	REAL(wp) :: pf ! sigma value
2402	!!----------------------------------------------------------------------
2403	!
2404	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2405	& - TANH( rn_thetb * rn_theta ) ) &
2406	& * ( COSH( rn_theta ) &
2407	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2408	& / ( 2._wp * SINH( rn_theta ) )
2409	!
2410	END FUNCTION fssig
2411
2412
2413	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2414	!!----------------------------------------------------------------------
2415	!! * ROUTINE fssig1 *
2416	!!
2417	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2418	!!
2419	!! ** Method : the function provides the non-dimensional position of
2420	!! T and W (i.e. between 0 and 1)
2421	!! T-points at integer values (between 1 and jpk)
2422	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2423	!!----------------------------------------------------------------------
2424	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2425	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2426	REAL(wp) :: pf1 ! sigma value
2427	!!----------------------------------------------------------------------
2428	!
2429	IF ( rn_theta == 0 ) then ! uniform sigma
2430	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2431	ELSE ! stretched sigma
2432	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2433	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2434	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2435	ENDIF
2436	!
2437	END FUNCTION fssig1
2438
2439
2440	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2441	!!----------------------------------------------------------------------
2442	!! * ROUTINE fgamma *
2443	!!
2444	!! ** Purpose : provide analytical function for the s-coordinate
2445	!!
2446	!! ** Method : the function provides the non-dimensional position of
2447	!! T and W (i.e. between 0 and 1)
2448	!! T-points at integer values (between 1 and jpk)
2449	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2450	!!
2451	!! This method allows the maintenance of fixed surface and or
2452	!! bottom cell resolutions (cf. geopotential coordinates)
2453	!! within an analytically derived stretched S-coordinate framework.
2454	!!
2455	!! Reference : Siddorn and Furner, in prep
2456	!!----------------------------------------------------------------------
2457	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2458	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2459	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2460	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2461	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2462	REAL(wp) :: za1,za2,za3 ! local variables
2463	REAL(wp) :: zn1,zn2 ! local variables
2464	REAL(wp) :: za,zb,zx ! local variables
2465	integer :: jk
2466	!!----------------------------------------------------------------------
2467	!
2468
2469	zn1 = 1./(jpk-1.)
2470	zn2 = 1. - zn1
2471
2472	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2473	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2474	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2475
2476	za = pzb - za3*(pzs-za1)-za2
2477	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2478	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2479	zx = 1.0_wp-za/2.0_wp-zb
2480
2481	DO jk = 1, jpk
2482	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2483	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2484	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2485	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2486	ENDDO
2487
2488	!
2489	END FUNCTION fgamma
2490
2491	!!======================================================================
2492	END MODULE domzgr

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

Download in other formats: