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

Context Navigation

source: branches/2013/dev_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4180

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

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

Download in other formats: