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/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5466

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

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

Download in other formats: