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

Last change on this file since 5074 was 5074, checked in by hliu, 9 years ago
prevent nla10 from being 0
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	DO ji = 1, jpi
1565	DO jj = 1, jpj
1566
1567	IF( bathy(ji,jj) > 0._wp) THEN
1568	DO jk = 1, mbathy(ji,jj)
1569	! check coordinate is monotonically increasing
1570	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1571	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1572	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1573	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
1574	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
1575	CALL ctl_stop( ctmp1 )
1576	ENDIF
1577	! and check it has never gone negative
1578	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1579	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1580	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1581	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1582	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1583	CALL ctl_stop( ctmp1 )
1584	ENDIF
1585	! and check it never exceeds the total depth
1586	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
1587	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1588	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1589	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1590	CALL ctl_stop( ctmp1 )
1591	ENDIF
1592	END DO
1593
1594	DO jk = 1, mbathy(ji,jj)-1
1595	! and check it never exceeds the total depth
1596	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
1597	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1598	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1599	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1600	CALL ctl_stop( ctmp1 )
1601	ENDIF
1602	END DO
1603
1604	ENDIF
1605
1606	END DO
1607	END DO
1608	!
1609	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1610	!
1611	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
1612	!
1613	END SUBROUTINE zgr_sco
1614
1615	!!======================================================================
1616	SUBROUTINE s_sh94()
1617
1618	!!----------------------------------------------------------------------
1619	!! * ROUTINE s_sh94 *
1620	!!
1621	!! ** Purpose : stretch the s-coordinate system
1622	!!
1623	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
1624	!! mixed S/sigma coordinate
1625	!!
1626	!! Reference : Song and Haidvogel 1994.
1627	!!----------------------------------------------------------------------
1628	!
1629	INTEGER :: ji, jj, jk ! dummy loop argument
1630	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1631	REAL(wp) :: ztmpu, ztmpv, ztmpf
1632	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1633	!
1634	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1635	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1636
1637	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1638	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1639
1640	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1641	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1642	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1643	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1644
1645	DO ji = 1, jpi
1646	DO jj = 1, jpj
1647
1648	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
1649	DO jk = 1, jpk
1650	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1651	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1652	END DO
1653	ELSE ! shallow water, uniform sigma
1654	DO jk = 1, jpk
1655	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
1656	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1657	END DO
1658	ENDIF
1659	!
1660	DO jk = 1, jpkm1
1661	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1662	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1663	END DO
1664	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
1665	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1666	!
1667	! Coefficients for vertical depth as the sum of e3w scale factors
1668	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
1669	DO jk = 2, jpk
1670	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1671	END DO
1672	!
1673	DO jk = 1, jpk
1674	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1675	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1676	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
1677	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
1678	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
1679	END DO
1680	!
1681	END DO ! for all jj's
1682	END DO ! for all ji's
1683
1684	DO ji = 1, jpim1
1685	DO jj = 1, jpjm1
1686	! extended for Wetting/Drying case
1687	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1688	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1689	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1690	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1691	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1692	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1693	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1694	DO jk = 1, jpk
1695	IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
1696	z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
1697	z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
1698	ELSE
1699	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1700	& / ztmpu
1701	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1702	& / ztmpu
1703	END IF
1704
1705	IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
1706	z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
1707	z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
1708	ELSE
1709	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1710	& / ztmpv
1711	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1712	& / ztmpv
1713	END IF
1714
1715	IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
1716	z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) &
1717	& + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
1718	ELSE
1719	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
1720	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
1721	& / ztmpf
1722	END IF
1723
1724	!
1725	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1726	! & / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1727	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1728	! & / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1729	!z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
1730	! & + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
1731	! & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1732	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1733	! & / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1734	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1735	! & / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1736	!
1737	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1738	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1739	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1740	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1741	!
1742	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1743	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1744	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1745	END DO
1746	END DO
1747	END DO
1748
1749	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1750	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1751
1752	END SUBROUTINE s_sh94
1753
1754	SUBROUTINE s_sf12
1755
1756	!!----------------------------------------------------------------------
1757	!! * ROUTINE s_sf12 *
1758	!!
1759	!! ** Purpose : stretch the s-coordinate system
1760	!!
1761	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
1762	!! mixed S/sigma/Z coordinate
1763	!!
1764	!! This method allows the maintenance of fixed surface and or
1765	!! bottom cell resolutions (cf. geopotential coordinates)
1766	!! within an analytically derived stretched S-coordinate framework.
1767	!!
1768	!!
1769	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1770	!!----------------------------------------------------------------------
1771	!
1772	INTEGER :: ji, jj, jk ! dummy loop argument
1773	REAL(wp) :: zsmth ! smoothing around critical depth
1774	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
1775	REAL(wp) :: ztmpu, ztmpv, ztmpf
1776	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1777	!
1778	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1779	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1780
1781	!
1782	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1783	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1784
1785	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1786	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1787	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1788	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1789
1790	DO ji = 1, jpi
1791	DO jj = 1, jpj
1792
1793	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1794
1795	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
1796	! could be changed by users but care must be taken to do so carefully
1797	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1798
1799	zzs = rn_zs / hbatt(ji,jj)
1800
1801	IF (rn_efold /= 0.0_wp) THEN
1802	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1803	ELSE
1804	zsmth = 1.0_wp
1805	ENDIF
1806
1807	DO jk = 1, jpk
1808	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1809	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1810	ENDDO
1811	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
1812	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
1813
1814	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1815
1816	DO jk = 1, jpk
1817	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1818	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1819	END DO
1820
1821	ELSE ! shallow water, z coordinates
1822
1823	DO jk = 1, jpk
1824	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1825	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1826	END DO
1827
1828	ENDIF
1829
1830	DO jk = 1, jpkm1
1831	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1832	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1833	END DO
1834	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
1835	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1836
1837	! Coefficients for vertical depth as the sum of e3w scale factors
1838	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
1839	DO jk = 2, jpk
1840	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1841	END DO
1842
1843	DO jk = 1, jpk
1844	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1845	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1846	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
1847	END DO
1848
1849	ENDDO ! for all jj's
1850	ENDDO ! for all ji's
1851
1852	DO ji=1,jpi-1
1853	DO jj=1,jpj-1
1854
1855	! extend to suit for Wetting/Drying case
1856	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1857	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1858	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1859	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1860	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1861	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1862	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1863	DO jk = 1, jpk
1864	IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN
1865	z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
1866	z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
1867	ELSE
1868	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1869	& / ztmpu
1870	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1871	& / ztmpu
1872	END IF
1873
1874	IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN
1875	z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
1876	z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
1877	ELSE
1878	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1879	& / ztmpv
1880	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1881	& / ztmpv
1882	END IF
1883
1884	IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN
1885	z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) &
1886	& + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
1887	ELSE
1888	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
1889	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
1890	& / ztmpf
1891	END IF
1892
1893	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
1894	! ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1895	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
1896	! ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1897	!z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
1898	! hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
1899	! ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1900	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
1901	! ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1902	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
1903	! ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1904
1905	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1906	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1907	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1908	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1909	!
1910	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1911	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1912	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1913	END DO
1914
1915	ENDDO
1916	ENDDO
1917	!
1918	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
1919	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
1920	CALL lbc_lnk(e3w_0 ,'T',1.)
1921	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
1922	!
1923	! ! =============
1924
1925	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1926	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1927
1928	END SUBROUTINE s_sf12
1929
1930	SUBROUTINE s_tanh()
1931
1932	!!----------------------------------------------------------------------
1933	!! * ROUTINE s_tanh*
1934	!!
1935	!! ** Purpose : stretch the s-coordinate system
1936	!!
1937	!! ** Method : s-coordinate stretch
1938	!!
1939	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1940	!!----------------------------------------------------------------------
1941
1942	INTEGER :: ji, jj, jk ! dummy loop argument
1943	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1944
1945	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
1946	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
1947
1948	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
1949	CALL wrk_alloc( jpk, z_esigt, z_esigw )
1950
1951	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
1952	z_esigt = 0._wp ; z_esigw = 0._wp
1953
1954	DO jk = 1, jpk
1955	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1956	z_gsigt(jk) = -fssig( REAL(jk,wp) )
1957	END DO
1958	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
1959	!
1960	! Coefficients for vertical scale factors at w-, t- levels
1961	!!gm bug : define it from analytical function, not like juste bellow....
1962	!!gm or betteroffer the 2 possibilities....
1963	DO jk = 1, jpkm1
1964	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
1965	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1966	END DO
1967	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
1968	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1969	!
1970	! Coefficients for vertical depth as the sum of e3w scale factors
1971	z_gsi3w(1) = 0.5_wp * z_esigw(1)
1972	DO jk = 2, jpk
1973	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
1974	END DO
1975	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1976	DO jk = 1, jpk
1977	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1978	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1979	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
1980	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
1981	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
1982	END DO
1983	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1984	DO jj = 1, jpj
1985	DO ji = 1, jpi
1986	DO jk = 1, jpk
1987	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1988	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1989	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1990	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1991	!
1992	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1993	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1994	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1995	END DO
1996	END DO
1997	END DO
1998
1999	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2000	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2001
2002	END SUBROUTINE s_tanh
2003
2004	FUNCTION fssig( pk ) RESULT( pf )
2005	!!----------------------------------------------------------------------
2006	!! * ROUTINE fssig *
2007	!!
2008	!! ** Purpose : provide the analytical function in s-coordinate
2009	!!
2010	!! ** Method : the function provide the non-dimensional position of
2011	!! T and W (i.e. between 0 and 1)
2012	!! T-points at integer values (between 1 and jpk)
2013	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2014	!!----------------------------------------------------------------------
2015	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2016	REAL(wp) :: pf ! sigma value
2017	!!----------------------------------------------------------------------
2018	!
2019	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2020	& - TANH( rn_thetb * rn_theta ) ) &
2021	& * ( COSH( rn_theta ) &
2022	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2023	& / ( 2._wp * SINH( rn_theta ) )
2024	!
2025	END FUNCTION fssig
2026
2027
2028	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2029	!!----------------------------------------------------------------------
2030	!! * ROUTINE fssig1 *
2031	!!
2032	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2033	!!
2034	!! ** Method : the function provides the non-dimensional position of
2035	!! T and W (i.e. between 0 and 1)
2036	!! T-points at integer values (between 1 and jpk)
2037	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2038	!!----------------------------------------------------------------------
2039	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2040	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2041	REAL(wp) :: pf1 ! sigma value
2042	!!----------------------------------------------------------------------
2043	!
2044	IF ( rn_theta == 0 ) then ! uniform sigma
2045	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2046	ELSE ! stretched sigma
2047	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2048	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2049	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2050	ENDIF
2051	!
2052	END FUNCTION fssig1
2053
2054
2055	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2056	!!----------------------------------------------------------------------
2057	!! * ROUTINE fgamma *
2058	!!
2059	!! ** Purpose : provide analytical function for the s-coordinate
2060	!!
2061	!! ** Method : the function provides the non-dimensional position of
2062	!! T and W (i.e. between 0 and 1)
2063	!! T-points at integer values (between 1 and jpk)
2064	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2065	!!
2066	!! This method allows the maintenance of fixed surface and or
2067	!! bottom cell resolutions (cf. geopotential coordinates)
2068	!! within an analytically derived stretched S-coordinate framework.
2069	!!
2070	!! Reference : Siddorn and Furner, in prep
2071	!!----------------------------------------------------------------------
2072	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2073	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2074	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2075	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2076	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2077	REAL(wp) :: za1,za2,za3 ! local variables
2078	REAL(wp) :: zn1,zn2 ! local variables
2079	REAL(wp) :: za,zb,zx ! local variables
2080	integer :: jk
2081	!!----------------------------------------------------------------------
2082	!
2083
2084	zn1 = 1./(jpk-1.)
2085	zn2 = 1. - zn1
2086
2087	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2088	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2089	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2090
2091	za = pzb - za3*(pzs-za1)-za2
2092	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2093	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2094	zx = 1.0_wp-za/2.0_wp-zb
2095
2096	DO jk = 1, jpk
2097	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2098	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2099	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2100	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2101	ENDDO
2102
2103	!
2104	END FUNCTION fgamma
2105
2106	!!======================================================================
2107	END MODULE domzgr

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

Download in other formats: