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/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3989

Last change on this file since 3989 was 3989, checked in by clevy, 11 years ago
Configuration setting/Step3 and doc, see ticket:#1074
Property svn:keywords set to `Id`
File size: 96.7 KB

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

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

Download in other formats: