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

Last change on this file since 4171 was 4171, checked in by rfurner, 11 years ago
merging changes from 3987:4141 of dev_r3987_UKMO4_OBS
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 = .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	!
1105	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1106
1107	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1108	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1109	!!----------------------------------------------------------------------
1110	!
1111	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1112	!
1113	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat )
1114	!
1115	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
1116	READ ( numnam, namzgr_sco )
1117
1118	IF(lwp) THEN ! control print
1119	WRITE(numout,*)
1120	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1121	WRITE(numout,*) '~~~~~~~~~~~'
1122	WRITE(numout,*) ' Namelist namzgr_sco'
1123	WRITE(numout,*) ' stretching coeffs '
1124	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1125	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1126	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1127	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1128	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1129	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1130	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1131	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1132	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1133	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1134	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1135	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1136	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1137	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1138	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1139	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1140	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1141	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1142	ENDIF
1143
1144	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1145	hifu(:,:) = rn_sbot_min
1146	hifv(:,:) = rn_sbot_min
1147	hiff(:,:) = rn_sbot_min
1148
1149	! ! set maximum ocean depth
1150	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1151
1152	DO jj = 1, jpj
1153	DO ji = 1, jpi
1154	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1155	END DO
1156	END DO
1157	! ! =============================
1158	! ! Define the envelop bathymetry (hbatt)
1159	! ! =============================
1160	! use r-value to create hybrid coordinates
1161	DO jj = 1, jpj
1162	DO ji = 1, jpi
1163	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1164	END DO
1165	END DO
1166	!
1167	! Smooth the bathymetry (if required)
1168	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1169	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1170	!
1171	jl = 0
1172	zrmax = 1._wp
1173	! ! ================ !
1174	DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
1175	! ! ================ !
1176	jl = jl + 1
1177	zrmax = 0._wp
1178	zmsk(:,:) = 0._wp
1179	DO jj = 1, nlcj
1180	DO ji = 1, nlci
1181	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1182	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1183	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1184	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1185	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1186	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp
1187	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp
1188	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp
1189	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp
1190	END DO
1191	END DO
1192	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1193	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1194	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )
1195	DO jj = 1, nlcj
1196	DO ji = 1, nlci
1197	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1198	END DO
1199	END DO
1200	!
1201	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1202	!
1203	DO jj = 1, nlcj
1204	DO ji = 1, nlci
1205	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1206	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1207	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1208	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1209	IF( zmsk(ji,jj) == 1._wp ) THEN
1210	ztmp(ji,jj) = ( &
1211	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1212	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) &
1213	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1214	& ) / ( &
1215	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1216	& + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) &
1217	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1218	& )
1219	ENDIF
1220	END DO
1221	END DO
1222	!
1223	DO jj = 1, nlcj
1224	DO ji = 1, nlci
1225	IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1226	END DO
1227	END DO
1228	!
1229	! Apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1230	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp )
1231	DO jj = 1, nlcj
1232	DO ji = 1, nlci
1233	IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj)
1234	END DO
1235	END DO
1236	! ! ================ !
1237	END DO ! End loop !
1238	! ! ================ !
1239	!
1240	! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions
1241	DO ji = nlci+1, jpi
1242	zenv(ji,1:nlcj) = zenv(nlci,1:nlcj)
1243	END DO
1244	!
1245	DO jj = nlcj+1, jpj
1246	zenv(:,jj) = zenv(:,nlcj)
1247	END DO
1248	!
1249	! Envelope bathymetry saved in hbatt
1250	hbatt(:,:) = zenv(:,:)
1251	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1252	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1253	DO jj = 1, jpj
1254	DO ji = 1, jpi
1255	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1256	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1257	END DO
1258	END DO
1259	ENDIF
1260	!
1261	IF(lwp) THEN ! Control print
1262	WRITE(numout,*)
1263	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1264	WRITE(numout,*)
1265	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1266	IF( nprint == 1 ) THEN
1267	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1268	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1269	ENDIF
1270	ENDIF
1271
1272	! ! ==============================
1273	! ! hbatu, hbatv, hbatf fields
1274	! ! ==============================
1275	IF(lwp) THEN
1276	WRITE(numout,*)
1277	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1278	ENDIF
1279	hbatu(:,:) = rn_sbot_min
1280	hbatv(:,:) = rn_sbot_min
1281	hbatf(:,:) = rn_sbot_min
1282	DO jj = 1, jpjm1
1283	DO ji = 1, jpim1 ! NO vector opt.
1284	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1285	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1286	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1287	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1288	END DO
1289	END DO
1290	!
1291	! Apply lateral boundary condition
1292	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1293	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
1294	DO jj = 1, jpj
1295	DO ji = 1, jpi
1296	IF( hbatu(ji,jj) == 0._wp ) THEN
1297	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1298	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1299	ENDIF
1300	END DO
1301	END DO
1302	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
1303	DO jj = 1, jpj
1304	DO ji = 1, jpi
1305	IF( hbatv(ji,jj) == 0._wp ) THEN
1306	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1307	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1308	ENDIF
1309	END DO
1310	END DO
1311	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
1312	DO jj = 1, jpj
1313	DO ji = 1, jpi
1314	IF( hbatf(ji,jj) == 0._wp ) THEN
1315	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1316	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1317	ENDIF
1318	END DO
1319	END DO
1320
1321	!!bug: key_helsinki a verifer
1322	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1323	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1324	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1325	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1326
1327	IF( nprint == 1 .AND. lwp ) THEN
1328	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1329	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1330	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1331	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1332	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1333	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1334	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1335	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1336	ENDIF
1337	!! helsinki
1338
1339	! ! =======================
1340	! ! s-ccordinate fields (gdep., e3.)
1341	! ! =======================
1342	!
1343	! non-dimensional "sigma" for model level depth at w- and t-levels
1344
1345
1346	!========================================================================
1347	! Song and Haidvogel 1994 (ln_s_sh94=T)
1348	! Siddorn and Furner 2012 (ln_sf12=T)
1349	! or tanh function (both false)
1350	!========================================================================
1351	IF ( ln_s_sh94 ) THEN
1352	CALL s_sh94()
1353	ELSE IF ( ln_s_sf12 ) THEN
1354	CALL s_sf12()
1355	ELSE
1356	CALL s_tanh()
1357	ENDIF
1358
1359	CALL lbc_lnk( e3t , 'T', 1._wp )
1360	CALL lbc_lnk( e3u , 'U', 1._wp )
1361	CALL lbc_lnk( e3v , 'V', 1._wp )
1362	CALL lbc_lnk( e3f , 'F', 1._wp )
1363	CALL lbc_lnk( e3w , 'W', 1._wp )
1364	CALL lbc_lnk( e3uw, 'U', 1._wp )
1365	CALL lbc_lnk( e3vw, 'V', 1._wp )
1366
1367	fsdepw(:,:,:) = gdepw (:,:,:)
1368	fsde3w(:,:,:) = gdep3w(:,:,:)
1369	!
1370	where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1._wp
1371	where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1._wp
1372	where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1._wp
1373	where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1._wp
1374	where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1._wp
1375	where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1._wp
1376	where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1._wp
1377
1378	#if defined key_agrif
1379	! Ensure meaningful vertical scale factors in ghost lines/columns
1380	IF( .NOT. Agrif_Root() ) THEN
1381	!
1382	IF((nbondi == -1).OR.(nbondi == 2)) THEN
1383	e3u(1,:,:) = e3u(2,:,:)
1384	ENDIF
1385	!
1386	IF((nbondi == 1).OR.(nbondi == 2)) THEN
1387	e3u(nlci-1,:,:) = e3u(nlci-2,:,:)
1388	ENDIF
1389	!
1390	IF((nbondj == -1).OR.(nbondj == 2)) THEN
1391	e3v(:,1,:) = e3v(:,2,:)
1392	ENDIF
1393	!
1394	IF((nbondj == 1).OR.(nbondj == 2)) THEN
1395	e3v(:,nlcj-1,:) = e3v(:,nlcj-2,:)
1396	ENDIF
1397	!
1398	ENDIF
1399	#endif
1400
1401	fsdept(:,:,:) = gdept (:,:,:)
1402	fsdepw(:,:,:) = gdepw (:,:,:)
1403	fsde3w(:,:,:) = gdep3w(:,:,:)
1404	fse3t (:,:,:) = e3t (:,:,:)
1405	fse3u (:,:,:) = e3u (:,:,:)
1406	fse3v (:,:,:) = e3v (:,:,:)
1407	fse3f (:,:,:) = e3f (:,:,:)
1408	fse3w (:,:,:) = e3w (:,:,:)
1409	fse3uw(:,:,:) = e3uw (:,:,:)
1410	fse3vw(:,:,:) = e3vw (:,:,:)
1411	!!
1412	! HYBRID :
1413	DO jj = 1, jpj
1414	DO ji = 1, jpi
1415	DO jk = 1, jpkm1
1416	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1417	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
1418	END DO
1419	END DO
1420	END DO
1421	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1422	& ' MAX ', MAXVAL( mbathy(:,:) )
1423
1424	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1425	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1426	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
1427	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
1428	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
1429	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
1430	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
1431	& ' w ', MINVAL( fse3w (:,:,:) )
1432
1433	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
1434	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
1435	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
1436	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
1437	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
1438	& ' w ', MAXVAL( fse3w (:,:,:) )
1439	ENDIF
1440	! END DO
1441	IF(lwp) THEN ! selected vertical profiles
1442	WRITE(numout,*)
1443	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1444	WRITE(numout,*) ' ~~~~~~ --------------------'
1445	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1446	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
1447	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1448	iip1 = MIN(20, jpiglo-1) ! for config with i smaller than 20 points
1449	ijp1 = MIN(20, jpjglo-1) ! for config with j smaller than 20 points
1450	DO jj = mj0(ijp1), mj1(ijp1)
1451	DO ji = mi0(iip1), mi1(iip1)
1452	WRITE(numout,*)
1453	WRITE(numout,*) ' domzgr: vertical coordinates : point (',iip1,',',ijp1,',k) bathy = ', &
1454	& bathy(ji,jj), hbatt(ji,jj)
1455	WRITE(numout,*) ' ~~~~~~ --------------------'
1456	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1457	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1458	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1459	END DO
1460	END DO
1461	iip1 = MIN( 74, jpiglo-1)
1462	ijp1 = MIN( 100, jpjglo-1)
1463	DO jj = mj0(ijp1), mj1(ijp1)
1464	DO ji = mi0(iip1), mi1(iip1)
1465	WRITE(numout,*)
1466	WRITE(numout,*) ' domzgr: vertical coordinates : point (',iip1,',',ijp1,',k) bathy = ', &
1467	& bathy(ji,jj), hbatt(ji,jj)
1468	WRITE(numout,*) ' ~~~~~~ --------------------'
1469	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1470	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1471	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1472	END DO
1473	END DO
1474	ENDIF
1475
1476	!================================================================================
1477	! check the coordinate makes sense
1478	!================================================================================
1479	DO ji = 1, jpi
1480	DO jj = 1, jpj
1481
1482	IF( hbatt(ji,jj) > 0._wp) THEN
1483	DO jk = 1, mbathy(ji,jj)
1484	! check coordinate is monotonically increasing
1485	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1486	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1487	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1488	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
1489	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
1490	CALL ctl_stop( ctmp1 )
1491	ENDIF
1492	! and check it has never gone negative
1493	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1494	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1495	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1496	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1497	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1498	CALL ctl_stop( ctmp1 )
1499	ENDIF
1500	! and check it never exceeds the total depth
1501	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
1502	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1503	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1504	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1505	CALL ctl_stop( ctmp1 )
1506	ENDIF
1507	END DO
1508
1509	DO jk = 1, mbathy(ji,jj)-1
1510	! and check it never exceeds the total depth
1511	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
1512	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1513	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1514	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1515	CALL ctl_stop( ctmp1 )
1516	ENDIF
1517	END DO
1518
1519	ENDIF
1520
1521	END DO
1522	END DO
1523	!
1524	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat )
1525	!
1526	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
1527	!
1528	END SUBROUTINE zgr_sco
1529
1530	!!======================================================================
1531	SUBROUTINE s_sh94()
1532
1533	!!----------------------------------------------------------------------
1534	!! * ROUTINE s_sh94 *
1535	!!
1536	!! ** Purpose : stretch the s-coordinate system
1537	!!
1538	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
1539	!! mixed S/sigma coordinate
1540	!!
1541	!! Reference : Song and Haidvogel 1994.
1542	!!----------------------------------------------------------------------
1543	!
1544	INTEGER :: ji, jj, jk ! dummy loop argument
1545	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1546	!
1547	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1548	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1549
1550	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1551	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1552
1553	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1554	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1555	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1556	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1557
1558	DO ji = 1, jpi
1559	DO jj = 1, jpj
1560
1561	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
1562	DO jk = 1, jpk
1563	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1564	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1565	END DO
1566	ELSE ! shallow water, uniform sigma
1567	DO jk = 1, jpk
1568	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
1569	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1570	END DO
1571	ENDIF
1572	!
1573	DO jk = 1, jpkm1
1574	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1575	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1576	END DO
1577	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
1578	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1579	!
1580	! Coefficients for vertical depth as the sum of e3w scale factors
1581	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
1582	DO jk = 2, jpk
1583	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1584	END DO
1585	!
1586	DO jk = 1, jpk
1587	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1588	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1589	gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
1590	gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
1591	gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
1592	END DO
1593	!
1594	END DO ! for all jj's
1595	END DO ! for all ji's
1596
1597	DO ji = 1, jpim1
1598	DO jj = 1, jpjm1
1599	DO jk = 1, jpk
1600	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1601	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1602	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1603	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1604	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
1605	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
1606	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1607	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1608	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1609	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1610	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1611	!
1612	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1613	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1614	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1615	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1616	!
1617	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1618	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1619	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1620	END DO
1621	END DO
1622	END DO
1623
1624	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1625	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1626
1627	END SUBROUTINE s_sh94
1628
1629	SUBROUTINE s_sf12
1630
1631	!!----------------------------------------------------------------------
1632	!! * ROUTINE s_sf12 *
1633	!!
1634	!! ** Purpose : stretch the s-coordinate system
1635	!!
1636	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
1637	!! mixed S/sigma/Z coordinate
1638	!!
1639	!! This method allows the maintenance of fixed surface and or
1640	!! bottom cell resolutions (cf. geopotential coordinates)
1641	!! within an analytically derived stretched S-coordinate framework.
1642	!!
1643	!!
1644	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1645	!!----------------------------------------------------------------------
1646	!
1647	INTEGER :: ji, jj, jk ! dummy loop argument
1648	REAL(wp) :: zsmth ! smoothing around critical depth
1649	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
1650	!
1651	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1652	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1653
1654	!
1655	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1656	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1657
1658	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1659	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1660	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1661	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1662
1663	DO ji = 1, jpi
1664	DO jj = 1, jpj
1665
1666	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1667
1668	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
1669	! could be changed by users but care must be taken to do so carefully
1670	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1671
1672	zzs = rn_zs / hbatt(ji,jj)
1673
1674	IF (rn_efold /= 0.0_wp) THEN
1675	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1676	ELSE
1677	zsmth = 1.0_wp
1678	ENDIF
1679
1680	DO jk = 1, jpk
1681	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1682	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1683	ENDDO
1684	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
1685	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
1686
1687	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1688
1689	DO jk = 1, jpk
1690	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1691	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1692	END DO
1693
1694	ELSE ! shallow water, z coordinates
1695
1696	DO jk = 1, jpk
1697	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1698	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1699	END DO
1700
1701	ENDIF
1702
1703	DO jk = 1, jpkm1
1704	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1705	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1706	END DO
1707	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
1708	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1709
1710	! Coefficients for vertical depth as the sum of e3w scale factors
1711	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
1712	DO jk = 2, jpk
1713	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1714	END DO
1715
1716	DO jk = 1, jpk
1717	gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1718	gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1719	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
1720	END DO
1721
1722	ENDDO ! for all jj's
1723	ENDDO ! for all ji's
1724
1725	DO ji=1,jpi-1
1726	DO jj=1,jpj-1
1727
1728	DO jk = 1, jpk
1729	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
1730	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1731	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
1732	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1733	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
1734	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
1735	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1736	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
1737	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1738	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
1739	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1740
1741	e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1742	e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1743	e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1744	e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1745	!
1746	e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1747	e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1748	e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1749	END DO
1750
1751	ENDDO
1752	ENDDO
1753	!
1754	! ! =============
1755
1756	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1757	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1758
1759	END SUBROUTINE s_sf12
1760
1761	SUBROUTINE s_tanh()
1762
1763	!!----------------------------------------------------------------------
1764	!! * ROUTINE s_tanh*
1765	!!
1766	!! ** Purpose : stretch the s-coordinate system
1767	!!
1768	!! ** Method : s-coordinate stretch
1769	!!
1770	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1771	!!----------------------------------------------------------------------
1772
1773	INTEGER :: ji, jj, jk ! dummy loop argument
1774	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1775
1776	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
1777	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
1778
1779	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
1780	CALL wrk_alloc( jpk, z_esigt, z_esigw )
1781
1782	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
1783	z_esigt = 0._wp ; z_esigw = 0._wp
1784
1785	DO jk = 1, jpk
1786	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1787	z_gsigt(jk) = -fssig( REAL(jk,wp) )
1788	END DO
1789	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
1790	!
1791	! Coefficients for vertical scale factors at w-, t- levels
1792	!!gm bug : define it from analytical function, not like juste bellow....
1793	!!gm or betteroffer the 2 possibilities....
1794	DO jk = 1, jpkm1
1795	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
1796	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1797	END DO
1798	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
1799	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1800	!
1801	! Coefficients for vertical depth as the sum of e3w scale factors
1802	z_gsi3w(1) = 0.5_wp * z_esigw(1)
1803	DO jk = 2, jpk
1804	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
1805	END DO
1806	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1807	DO jk = 1, jpk
1808	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1809	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1810	gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
1811	gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
1812	gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
1813	END DO
1814	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1815	DO jj = 1, jpj
1816	DO ji = 1, jpi
1817	DO jk = 1, jpk
1818	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1819	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1820	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1821	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1822	!
1823	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1824	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1825	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1826	END DO
1827	END DO
1828	END DO
1829
1830	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
1831	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
1832
1833	END SUBROUTINE s_tanh
1834
1835	FUNCTION fssig( pk ) RESULT( pf )
1836	!!----------------------------------------------------------------------
1837	!! * ROUTINE fssig *
1838	!!
1839	!! ** Purpose : provide the analytical function in s-coordinate
1840	!!
1841	!! ** Method : the function provide the non-dimensional position of
1842	!! T and W (i.e. between 0 and 1)
1843	!! T-points at integer values (between 1 and jpk)
1844	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1845	!!----------------------------------------------------------------------
1846	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
1847	REAL(wp) :: pf ! sigma value
1848	!!----------------------------------------------------------------------
1849	!
1850	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1,wp) + rn_thetb ) ) &
1851	& - TANH( rn_thetb * rn_theta ) ) &
1852	& * ( COSH( rn_theta ) &
1853	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
1854	& / ( 2._wp * SINH( rn_theta ) )
1855	!
1856	END FUNCTION fssig
1857
1858
1859	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1860	!!----------------------------------------------------------------------
1861	!! * ROUTINE fssig1 *
1862	!!
1863	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
1864	!!
1865	!! ** Method : the function provides the non-dimensional position of
1866	!! T and W (i.e. between 0 and 1)
1867	!! T-points at integer values (between 1 and jpk)
1868	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1869	!!----------------------------------------------------------------------
1870	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
1871	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
1872	REAL(wp) :: pf1 ! sigma value
1873	!!----------------------------------------------------------------------
1874	!
1875	IF ( rn_theta == 0 ) then ! uniform sigma
1876	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1,wp )
1877	ELSE ! stretched sigma
1878	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1,wp)) ) ) / SINH( rn_theta ) &
1879	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1,wp)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
1880	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
1881	ENDIF
1882	!
1883	END FUNCTION fssig1
1884
1885
1886	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
1887	!!----------------------------------------------------------------------
1888	!! * ROUTINE fgamma *
1889	!!
1890	!! ** Purpose : provide analytical function for the s-coordinate
1891	!!
1892	!! ** Method : the function provides the non-dimensional position of
1893	!! T and W (i.e. between 0 and 1)
1894	!! T-points at integer values (between 1 and jpk)
1895	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1896	!!
1897	!! This method allows the maintenance of fixed surface and or
1898	!! bottom cell resolutions (cf. geopotential coordinates)
1899	!! within an analytically derived stretched S-coordinate framework.
1900	!!
1901	!! Reference : Siddorn and Furner, in prep
1902	!!----------------------------------------------------------------------
1903	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
1904	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
1905	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
1906	REAL(wp), INTENT(in ) :: pzs ! surface box depth
1907	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
1908	REAL(wp) :: za1,za2,za3 ! local variables
1909	REAL(wp) :: zn1,zn2 ! local variables
1910	REAL(wp) :: za,zb,zx ! local variables
1911	integer :: jk
1912	!!----------------------------------------------------------------------
1913	!
1914
1915	zn1 = 1./(jpk-1.)
1916	zn2 = 1. - zn1
1917
1918	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
1919	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
1920	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
1921
1922	za = pzb - za3*(pzs-za1)-za2
1923	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
1924	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
1925	zx = 1.0_wp-za/2.0_wp-zb
1926
1927	DO jk = 1, jpk
1928	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
1929	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
1930	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
1931	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
1932	ENDDO
1933
1934	!
1935	END FUNCTION fgamma
1936
1937	!!======================================================================
1938	END MODULE domzgr

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

Download in other formats: