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

Context Navigation

source: trunk/NEMO/OPA_SRC/DOM/domzgr.F90 @ 1273

Last change on this file since 1273 was 1273, checked in by ctlod, 15 years ago
update Gibrlatar, Bab El Mandeb and Sound straits in both full & partial steps bathymetry files such as closed seas, see ticket: #305
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 68.2 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	!!----------------------------------------------------------------------
16
17	!!----------------------------------------------------------------------
18	!! dom_zgr : defined the ocean vertical coordinate system
19	!! zgr_bat : bathymetry fields (levels and meters)
20	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
21	!! zgr_bat_ctl : check the bathymetry files
22	!! zgr_z : reference z-coordinate
23	!! zgr_zco : z-coordinate
24	!! zgr_zps : z-coordinate with partial steps
25	!! zgr_sco : s-coordinate
26	!! fssig : sigma coordinate non-dimensional function
27	!! dfssig : derivative of the sigma coordinate function !!gm (currently missing!)
28	!!---------------------------------------------------------------------
29	USE oce ! ocean dynamics and tracers
30	USE dom_oce ! ocean space and time domain
31	USE in_out_manager ! I/O manager
32	USE lib_mpp ! distributed memory computing library
33	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
34	USE closea ! closed seas
35	USE solisl ! solver - island in rigid-lid
36	USE c1d ! 1D configuration
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC dom_zgr ! called by dom_init.F90
42
43	!!gm DOCTOR name for the namelist parameter...
44	! !!! Namelist nam_zgr_sco
45	REAL(wp) :: sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m)
46	REAL(wp) :: sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
47	REAL(wp) :: theta = 6.0 ! surface control parameter (0<=theta<=20)
48	REAL(wp) :: thetb = 0.75 ! bottom control parameter (0<=thetb<= 1)
49	REAL(wp) :: r_max = 0.15 ! maximum cut-off r-value allowed (0<r_max<1)
50
51	!! * Substitutions
52	# include "domzgr_substitute.h90"
53	# include "vectopt_loop_substitute.h90"
54	!!----------------------------------------------------------------------
55	!! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
56	!! $Id$
57	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
58	!!----------------------------------------------------------------------
59
60	CONTAINS
61
62	SUBROUTINE dom_zgr
63	!!----------------------------------------------------------------------
64	!! * ROUTINE dom_zgr *
65	!!
66	!! ** Purpose : set the depth of model levels and the resulting
67	!! vertical scale factors.
68	!!
69	!! ** Method : - reference 1D vertical coordinate (gdep._0, e3._0)
70	!! - read/set ocean depth and ocean levels (bathy, mbathy)
71	!! - vertical coordinate (gdep., e3.) depending on the
72	!! coordinate chosen :
73	!! ln_zco=T z-coordinate (forced if lk_zco)
74	!! ln_zps=T z-coordinate with partial steps
75	!! ln_zco=T s-coordinate
76	!!
77	!! ** Action : define gdep., e3., mbathy and bathy
78	!!----------------------------------------------------------------------
79	INTEGER :: ioptio = 0 ! temporary integer
80	!!
81	NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
82	!!----------------------------------------------------------------------
83
84	REWIND ( numnam ) ! Read Namelist nam_zgr : vertical coordinate'
85	READ ( numnam, nam_zgr )
86
87	IF(lwp) THEN ! Control print
88	WRITE(numout,*)
89	WRITE(numout,*) 'dom_zgr : vertical coordinate'
90	WRITE(numout,*) '~~~~~~~'
91	WRITE(numout,*) ' Namelist nam_zgr : set vertical coordinate'
92	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
93	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
94	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
95	ENDIF
96
97	ioptio = 0 ! Check Vertical coordinate options
98	IF( ln_zco ) ioptio = ioptio + 1
99	IF( ln_zps ) ioptio = ioptio + 1
100	IF( ln_sco ) ioptio = ioptio + 1
101	IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
102	IF( lk_zco ) THEN
103	IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement'
104	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
105	ENDIF
106
107	! Build the vertical coordinate system
108	! ------------------------------------
109	CALL zgr_z ! Reference z-coordinate system (always called)
110	CALL zgr_bat ! Bathymetry fields (levels and meters)
111	IF( ln_zco ) CALL zgr_zco ! z-coordinate
112	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
113	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
114	!
115	END SUBROUTINE dom_zgr
116
117
118	SUBROUTINE zgr_z
119	!!----------------------------------------------------------------------
120	!! * ROUTINE zgr_z *
121	!!
122	!! ** Purpose : set the depth of model levels and the resulting
123	!! vertical scale factors.
124	!!
125	!! ** Method : z-coordinate system (use in all type of coordinate)
126	!! The depth of model levels is defined from an analytical
127	!! function the derivative of which gives the scale factors.
128	!! both depth and scale factors only depend on k (1d arrays).
129	!! w-level: gdepw_0 = fsdep(k)
130	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
131	!! t-level: gdept_0 = fsdep(k+0.5)
132	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
133	!!
134	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
135	!! - e3t_0 , e3w_0 : scale factors at T- and W-levels (m)
136	!!
137	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
138	!!----------------------------------------------------------------------
139	INTEGER :: jk ! dummy loop indices
140	REAL(wp) :: zt, zw ! temporary scalars
141	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
142	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
143	!!----------------------------------------------------------------------
144
145	! Set variables from parameters
146	! ------------------------------
147	zkth = ppkth ; zacr = ppacr
148	zdzmin = ppdzmin ; zhmax = pphmax
149
150	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
151	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
152	!
153	IF( ppa1 == pp_to_be_computed .AND. &
154	& ppa0 == pp_to_be_computed .AND. &
155	& ppsur == pp_to_be_computed ) THEN
156	!
157	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
158	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
159	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
160	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
161	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
162	ELSE
163	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
164	ENDIF
165
166	IF(lwp) THEN ! Parameter print
167	WRITE(numout,*)
168	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
169	WRITE(numout,*) ' ~~~~~~~'
170	IF( ppkth == 0. ) THEN
171	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
172	WRITE(numout,*) ' Total depth :', zhmax
173	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
174	ELSE
175	IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
176	WRITE(numout,*) ' zsur, za0, za1 computed from '
177	WRITE(numout,*) ' zdzmin = ', zdzmin
178	WRITE(numout,*) ' zhmax = ', zhmax
179	ENDIF
180	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
181	WRITE(numout,*) ' zsur = ', zsur
182	WRITE(numout,*) ' za0 = ', za0
183	WRITE(numout,*) ' za1 = ', za1
184	WRITE(numout,*) ' zkth = ', zkth
185	WRITE(numout,*) ' zacr = ', zacr
186	ENDIF
187	ENDIF
188
189
190	! Reference z-coordinate (depth - scale factor at T- and W-points)
191	! ======================
192	IF( ppkth == 0.e0 ) THEN ! uniform vertical grid
193	za1 = zhmax / FLOAT(jpk-1)
194	DO jk = 1, jpk
195	zw = FLOAT( jk )
196	zt = FLOAT( jk ) + 0.5
197	gdepw_0(jk) = ( zw - 1 ) * za1
198	gdept_0(jk) = ( zt - 1 ) * za1
199	e3w_0 (jk) = za1
200	e3t_0 (jk) = za1
201	END DO
202	ELSE ! Madec & Imbard 1996 function
203	DO jk = 1, jpk
204	zw = FLOAT( jk )
205	zt = FLOAT( jk ) + 0.5
206	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
207	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
208	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
209	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
210	END DO
211	gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero
212	ENDIF
213
214	IF(lwp) THEN ! control print
215	WRITE(numout,*)
216	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
217	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
218	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
219	ENDIF
220	DO jk = 1, jpk ! control positivity
221	IF( e3w_0 (jk) <= 0.e0 .OR. e3t_0 (jk) <= 0.e0 ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
222	IF( gdepw_0(jk) < 0.e0 .OR. gdept_0(jk) < 0.e0 ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
223	END DO
224	!
225	END SUBROUTINE zgr_z
226
227
228	SUBROUTINE zgr_bat
229	!!----------------------------------------------------------------------
230	!! * ROUTINE zgr_bat *
231	!!
232	!! ** Purpose : set bathymetry both in levels and meters
233	!!
234	!! ** Method : read or define mbathy and bathy arrays
235	!! * level bathymetry:
236	!! The ocean basin geometry is given by a two-dimensional array,
237	!! mbathy, which is defined as follow :
238	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
239	!! at t-point (ji,jj).
240	!! = 0 over the continental t-point.
241	!! = -n over the nth island t-point.
242	!! The array mbathy is checked to verified its consistency with
243	!! model option. in particular:
244	!! mbathy must have at least 1 land grid-points (mbathy<=0)
245	!! along closed boundary.
246	!! mbathy must be cyclic IF jperio=1.
247	!! mbathy must be lower or equal to jpk-1.
248	!! isolated ocean grid points are suppressed from mbathy
249	!! since they are only connected to remaining
250	!! ocean through vertical diffusion.
251	!! ntopo=-1 : rectangular channel or bassin with a bump
252	!! ntopo= 0 : flat rectangular channel or basin
253	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
254	!! bathy is read in 'bathy_meter.nc' NetCDF file
255	!! C A U T I O N : mbathy will be modified during the initializa-
256	!! tion phase to become the number of non-zero w-levels of a water
257	!! column, with a minimum value of 1.
258	!!
259	!! ** Action : - mbathy: level bathymetry (in level index)
260	!! - bathy : meter bathymetry (in meters)
261	!!----------------------------------------------------------------------
262	USE iom
263	!!
264	INTEGER :: ji, jj, jl, jk ! dummy loop indices
265	INTEGER :: inum ! temporary logical unit
266	INTEGER :: ii0, ii1, ij0, ij1 ! indices
267	INTEGER :: ii_bump, ij_bump, ih ! bump center position
268	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
269	REAL(wp) :: zi , zj , zh ! temporary scalars
270	INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data
271	REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data
272	!!----------------------------------------------------------------------
273
274	IF(lwp) WRITE(numout,*)
275	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
276	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
277
278	! ! ================== !
279	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
280	! ! ================== !
281	! ! global domain level and meter bathymetry (idta,zdta)
282	!
283	IF( ntopo == 0 ) THEN ! flat basin
284	IF(lwp) WRITE(numout,*)
285	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
286	idta(:,:) = jpkm1 ! before last level
287	zdta(:,:) = gdepw_0(jpk) ! last w-point depth
288	h_oce = gdepw_0(jpk)
289	ELSE ! bump centered in the basin
290	IF(lwp) WRITE(numout,*)
291	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
292	ii_bump = jpidta / 2 ! i-index of the bump center
293	ij_bump = jpjdta / 2 ! j-index of the bump center
294	r_bump = 50000.e0 ! bump radius (meters)
295	h_bump = 2700.e0 ! bump height (meters)
296	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
297	IF(lwp) WRITE(numout,*) ' bump characteristics: '
298	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
299	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
300	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
301	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
302	!
303	DO jj = 1, jpjdta ! zdta :
304	DO ji = 1, jpidta
305	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
306	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
307	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
308	END DO
309	END DO
310	! ! idta :
311	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
312	idta(:,:) = jpkm1
313	ELSE ! z-coordinate (zco or zps): step-like topography
314	idta(:,:) = jpkm1
315	DO jk = 1, jpkm1
316	DO jj = 1, jpjdta
317	DO ji = 1, jpidta
318	IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk
319	END DO
320	END DO
321	END DO
322	ENDIF
323	ENDIF
324	! ! set GLOBAL boundary conditions
325	! ! Caution : idta on the global domain: use of jperio, not nperio
326	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
327	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0
328	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
329	ELSEIF( jperio == 2 ) THEN
330	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
331	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
332	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0
333	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0
334	ELSE
335	ih = 0 ; zh = 0.e0
336	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
337	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
338	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
339	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
340	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
341	ENDIF
342
343	! ! local domain level and meter bathymetries (mbathy,bathy)
344	mbathy(:,:) = 0 ! set to zero extra halo points
345	bathy (:,:) = 0.e0 ! (require for mpp case)
346	DO jj = 1, nlcj ! interior values
347	DO ji = 1, nlci
348	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
349	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
350	END DO
351	END DO
352	!
353	! ! ================ !
354	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
355	! ! ================ !
356	!
357	IF( ln_zco ) THEN ! zco : read level bathymetry
358	CALL iom_open( 'bathy_level.nc', inum )
359	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
360	CALL iom_close (inum)
361	mbathy(:,:) = INT( bathy(:,:) )
362	! ! =====================
363	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
364	! ! =====================
365	IF( n_cla == 0 ) THEN
366	!
367	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
368	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
369	DO ji = mi0(ii0), mi1(ii1)
370	DO jj = mj0(ij0), mj1(ij1)
371	mbathy(ji,jj) = 15
372	END DO
373	END DO
374	IF(lwp) WRITE(numout,*)
375	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
376	!
377	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
378	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
379	DO ji = mi0(ii0), mi1(ii1)
380	DO jj = mj0(ij0), mj1(ij1)
381	mbathy(ji,jj) = 12
382	END DO
383	END DO
384	IF(lwp) WRITE(numout,*)
385	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
386	!
387	ENDIF
388	!
389	ENDIF
390	!
391	ENDIF
392	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
393	CALL iom_open( 'bathy_meter.nc', inum )
394	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
395	CALL iom_close (inum)
396	! ! =====================
397	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
398	! ! =====================
399	IF( n_cla == 0 ) THEN
400	!
401	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
402	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
403	DO ji = mi0(ii0), mi1(ii1)
404	DO jj = mj0(ij0), mj1(ij1)
405	bathy(ji,jj) = 284.e0
406	END DO
407	END DO
408	IF(lwp) WRITE(numout,*)
409	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
410	!
411	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
412	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
413	DO ji = mi0(ii0), mi1(ii1)
414	DO jj = mj0(ij0), mj1(ij1)
415	bathy(ji,jj) = 137.e0
416	END DO
417	END DO
418	IF(lwp) WRITE(numout,*)
419	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
420	!
421	ENDIF
422	!
423	ENDIF
424
425	ENDIF
426	! ! =============== !
427	ELSE ! error !
428	! ! =============== !
429	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
430	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
431	ENDIF
432	!
433	! ! =========================== !
434	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
435	DO jl = 1, jpncs ! =========================== !
436	DO jj = ncsj1(jl), ncsj2(jl)
437	DO ji = ncsi1(jl), ncsi2(jl)
438	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
439	bathy (ji,jj) = 0.e0
440	END DO
441	END DO
442	END DO
443	ENDIF
444	#if defined key_orca_lev10
445	! ! ================= !
446	mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !
447	! ! ================= !
448	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
449	#endif
450	! ! =============== !
451	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
452	! ! =============== !
453
454	! ! =================== !
455	IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check !
456	! ! =================== !
457	END SUBROUTINE zgr_bat
458
459
460	SUBROUTINE zgr_bat_zoom
461	!!----------------------------------------------------------------------
462	!! * ROUTINE zgr_bat_zoom *
463	!!
464	!! ** Purpose : - Close zoom domain boundary if necessary
465	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
466	!!
467	!! ** Method :
468	!!
469	!! ** Action : - update mbathy: level bathymetry (in level index)
470	!!----------------------------------------------------------------------
471	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
472	!!----------------------------------------------------------------------
473	!
474	IF(lwp) WRITE(numout,*)
475	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
476	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
477	!
478	! Zoom domain
479	! ===========
480	!
481	! Forced closed boundary if required
482	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
483	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
484	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
485	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
486	!
487	! Configuration specific domain modifications
488	! (here, ORCA arctic configuration: suppress Med Sea)
489	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
490	SELECT CASE ( jp_cfg )
491	! ! =======================
492	CASE ( 2 ) ! ORCA_R2 configuration
493	! ! =======================
494	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
495	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
496	ij0 = 98 ; ij1 = 110
497	! ! =======================
498	CASE ( 05 ) ! ORCA_R05 configuration
499	! ! =======================
500	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
501	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
502	ij0 = 314 ; ij1 = 370
503	END SELECT
504	!
505	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
506	!
507	ENDIF
508	!
509	END SUBROUTINE zgr_bat_zoom
510
511
512	SUBROUTINE zgr_bat_ctl
513	!!----------------------------------------------------------------------
514	!! * ROUTINE zgr_bat_ctl *
515	!!
516	!! ** Purpose : check the bathymetry in levels
517	!!
518	!! ** Method : The array mbathy is checked to verified its consistency
519	!! with the model options. in particular:
520	!! mbathy must have at least 1 land grid-points (mbathy<=0)
521	!! along closed boundary.
522	!! mbathy must be cyclic IF jperio=1.
523	!! mbathy must be lower or equal to jpk-1.
524	!! isolated ocean grid points are suppressed from mbathy
525	!! since they are only connected to remaining
526	!! ocean through vertical diffusion.
527	!! C A U T I O N : mbathy will be modified during the initializa-
528	!! tion phase to become the number of non-zero w-levels of a water
529	!! column, with a minimum value of 1.
530	!!
531	!! ** Action : - update mbathy: level bathymetry (in level index)
532	!! - update bathy : meter bathymetry (in meters)
533	!!----------------------------------------------------------------------
534	INTEGER :: ji, jj, jl ! dummy loop indices
535	INTEGER :: icompt, ibtest, ikmax ! temporary integers
536	REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace
537	!!----------------------------------------------------------------------
538
539	IF(lwp) WRITE(numout,*)
540	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
541	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
542
543	! ! Suppress isolated ocean grid points
544	IF(lwp) WRITE(numout,*)
545	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
546	IF(lwp) WRITE(numout,*)' -----------------------------------'
547	icompt = 0
548	DO jl = 1, 2
549	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
550	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
551	mbathy(jpi,:) = mbathy( 2 ,:)
552	ENDIF
553	DO jj = 2, jpjm1
554	DO ji = 2, jpim1
555	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
556	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
557	IF( ibtest < mbathy(ji,jj) ) THEN
558	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
559	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
560	mbathy(ji,jj) = ibtest
561	icompt = icompt + 1
562	ENDIF
563	END DO
564	END DO
565	END DO
566	IF( icompt == 0 ) THEN
567	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
568	ELSE
569	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
570	ENDIF
571	IF( lk_mpp ) THEN
572	zbathy(:,:) = FLOAT( mbathy(:,:) )
573	CALL lbc_lnk( zbathy, 'T', 1. )
574	mbathy(:,:) = INT( zbathy(:,:) )
575	ENDIF
576
577	! ! East-west cyclic boundary conditions
578	IF( nperio == 0 ) THEN
579	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
580	IF( lk_mpp ) THEN
581	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
582	IF( jperio /= 1 ) mbathy(1,:) = 0
583	ENDIF
584	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
585	IF( jperio /= 1 ) mbathy(nlci,:) = 0
586	ENDIF
587	ELSE
588	IF( ln_zco .OR. ln_zps ) THEN
589	mbathy( 1 ,:) = 0
590	mbathy(jpi,:) = 0
591	ELSE
592	mbathy( 1 ,:) = jpkm1
593	mbathy(jpi,:) = jpkm1
594	ENDIF
595	ENDIF
596	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
597	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
598	mbathy( 1 ,:) = mbathy(jpim1,:)
599	mbathy(jpi,:) = mbathy( 2 ,:)
600	ELSEIF( nperio == 2 ) THEN
601	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
602	ELSE
603	IF(lwp) WRITE(numout,*) ' e r r o r'
604	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
605	! STOP 'dom_mba'
606	ENDIF
607
608	! Set to zero mbathy over islands if necessary (lk_isl=F)
609	IF( .NOT. lk_isl ) THEN ! No island
610	IF(lwp) WRITE(numout,*)
611	IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'
612	IF(lwp) WRITE(numout,*) ' ----------------------------'
613	!
614	mbathy(:,:) = MAX( 0, mbathy(:,:) )
615	!
616	! Boundary condition on mbathy
617	IF( .NOT.lk_mpp ) THEN
618	!!gm !!bug ??? think about it !
619	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
620	zbathy(:,:) = FLOAT( mbathy(:,:) )
621	CALL lbc_lnk( zbathy, 'T', 1. )
622	mbathy(:,:) = INT( zbathy(:,:) )
623	ENDIF
624	ENDIF
625
626	! Number of ocean level inferior or equal to jpkm1
627	ikmax = 0
628	DO jj = 1, jpj
629	DO ji = 1, jpi
630	ikmax = MAX( ikmax, mbathy(ji,jj) )
631	END DO
632	END DO
633	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
634	IF( ikmax > jpkm1 ) THEN
635	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
636	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
637	ELSE IF( ikmax < jpkm1 ) THEN
638	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
639	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
640	ENDIF
641
642	! control print
643	IF( lwp .AND. nprint == 1 ) THEN
644	WRITE(numout,*)
645	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
646	WRITE(numout,*) ' ------------------'
647	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
648	WRITE(numout,*)
649	ENDIF
650	!
651	END SUBROUTINE zgr_bat_ctl
652
653
654	SUBROUTINE zgr_zco
655	!!----------------------------------------------------------------------
656	!! * ROUTINE zgr_zco *
657	!!
658	!! ** Purpose : define the z-coordinate system
659	!!
660	!! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T
661	!!----------------------------------------------------------------------
662	INTEGER :: jk
663	!!----------------------------------------------------------------------
664	!
665	IF( .NOT.lk_zco ) THEN
666	DO jk = 1, jpk
667	fsdept(:,:,jk) = gdept_0(jk)
668	fsdepw(:,:,jk) = gdepw_0(jk)
669	fsde3w(:,:,jk) = gdepw_0(jk)
670	fse3t (:,:,jk) = e3t_0(jk)
671	fse3u (:,:,jk) = e3t_0(jk)
672	fse3v (:,:,jk) = e3t_0(jk)
673	fse3f (:,:,jk) = e3t_0(jk)
674	fse3w (:,:,jk) = e3w_0(jk)
675	fse3uw(:,:,jk) = e3w_0(jk)
676	fse3vw(:,:,jk) = e3w_0(jk)
677	END DO
678	ENDIF
679	!
680	END SUBROUTINE zgr_zco
681
682	#if defined key_zco
683	!!----------------------------------------------------------------------
684	!! 'key_zco' : "pure" zco (gdep & e3 are 1D arrays)
685	!!----------------------------------------------------------------------
686	SUBROUTINE zgr_zps ! Empty routine
687	END SUBROUTINE zgr_zps
688	SUBROUTINE zgr_sco ! Empty routine
689	END SUBROUTINE zgr_sco
690
691	#else
692	!!----------------------------------------------------------------------
693	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
694	!!----------------------------------------------------------------------
695
696	SUBROUTINE zgr_zps
697	!!----------------------------------------------------------------------
698	!! * ROUTINE zgr_zps *
699	!!
700	!! ** Purpose : the depth and vertical scale factor in partial step
701	!! z-coordinate case
702	!!
703	!! ** Method : Partial steps : computes the 3D vertical scale factors
704	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
705	!! a partial step representation of bottom topography.
706	!!
707	!! The reference depth of model levels is defined from an analytical
708	!! function the derivative of which gives the reference vertical
709	!! scale factors.
710	!! From depth and scale factors reference, we compute there new value
711	!! with partial steps on 3d arrays ( i, j, k ).
712	!!
713	!! w-level: gdepw(i,j,k) = fsdep(k)
714	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
715	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
716	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
717	!!
718	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
719	!! we find the mbathy index of the depth at each grid point.
720	!! This leads us to three cases:
721	!!
722	!! - bathy = 0 => mbathy = 0
723	!! - 1 < mbathy < jpkm1
724	!! - bathy > gdepw(jpk) => mbathy = jpkm1
725	!!
726	!! Then, for each case, we find the new depth at t- and w- levels
727	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
728	!! and f-points.
729	!!
730	!! This routine is given as an example, it must be modified
731	!! following the user s desiderata. nevertheless, the output as
732	!! well as the way to compute the model levels and scale factors
733	!! must be respected in order to insure second order accuracy
734	!! schemes.
735	!!
736	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
737	!! - - - - - - - gdept, gdepw and e3. are positives
738	!!
739	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
740	!!----------------------------------------------------------------------
741	INTEGER :: ji, jj, jk ! dummy loop indices
742	INTEGER :: ik, it ! temporary integers
743	LOGICAL :: ll_print ! Allow control print for debugging
744	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
745	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
746	REAL(wp) :: zmax, zmin ! Maximum and minimum depth
747	REAL(wp) :: zdiff ! temporary scalar
748	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace
749	!!---------------------------------------------------------------------
750
751	IF(lwp) WRITE(numout,*)
752	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
753	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
754	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
755
756	ll_print=.FALSE. ! Local variable for debugging
757	!! ll_print=.TRUE.
758
759	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
760	WRITE(numout,*)
761	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
762	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
763	ENDIF
764
765
766	! bathymetry in level (from bathy_meter)
767	! ===================
768	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
769	zmin = gdepw_0(4) ! minimum depth = 3 levels
770
771	mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available
772
773	! ! storage of land and island's number (zera and negative values) in mbathy
774	WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )
775
776	! bounded value of bathy
777	!!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0
778	DO jj = 1, jpj
779	DO ji= 1, jpi
780	IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0
781	ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )
782	ENDIF
783	END DO
784	END DO
785
786	! Compute mbathy for ocean points (i.e. the number of ocean levels)
787	! find the number of ocean levels such that the last level thickness
788	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
789	! e3t_0 is the reference level thickness
790	DO jk = jpkm1, 1, -1
791	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
792	DO jj = 1, jpj
793	DO ji = 1, jpi
794	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
795	END DO
796	END DO
797	END DO
798
799	! Scale factors and depth at T- and W-points
800	DO jk = 1, jpk ! intitialization to the reference z-coordinate
801	gdept(:,:,jk) = gdept_0(jk)
802	gdepw(:,:,jk) = gdepw_0(jk)
803	e3t (:,:,jk) = e3t_0 (jk)
804	e3w (:,:,jk) = e3w_0 (jk)
805	END DO
806	hdept(:,:) = gdept(:,:,2 )
807	hdepw(:,:) = gdepw(:,:,3 )
808	!
809	DO jj = 1, jpj
810	DO ji = 1, jpi
811	ik = mbathy(ji,jj)
812	IF( ik > 0 ) THEN ! ocean point only
813	! max ocean level case
814	IF( ik == jpkm1 ) THEN
815	zdepwp = bathy(ji,jj)
816	ze3tp = bathy(ji,jj) - gdepw_0(ik)
817	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
818	e3t(ji,jj,ik ) = ze3tp
819	e3t(ji,jj,ik+1) = ze3tp
820	e3w(ji,jj,ik ) = ze3wp
821	e3w(ji,jj,ik+1) = ze3tp
822	gdepw(ji,jj,ik+1) = zdepwp
823	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
824	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
825	!
826	ELSE ! standard case
827	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
828	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
829	ENDIF
830	!gm Bug? check the gdepw_0
831	! ... on ik
832	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
833	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
834	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
835	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
836	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
837	e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
838	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
839	! ... on ik+1
840	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
841	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
842	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
843	ENDIF
844	ENDIF
845	END DO
846	END DO
847	!
848	it = 0
849	DO jj = 1, jpj
850	DO ji = 1, jpi
851	ik = mbathy(ji,jj)
852	IF( ik > 0 ) THEN ! ocean point only
853	hdept(ji,jj) = gdept(ji,jj,ik )
854	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
855	e3tp (ji,jj) = e3t(ji,jj,ik )
856	e3wp (ji,jj) = e3w(ji,jj,ik )
857	! test
858	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
859	IF( zdiff <= 0. .AND. lwp ) THEN
860	it = it + 1
861	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
862	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
863	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
864	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
865	ENDIF
866	ENDIF
867	END DO
868	END DO
869
870	! Scale factors and depth at U-, V-, UW and VW-points
871	DO jk = 1, jpk ! initialisation to z-scale factors
872	e3u (:,:,jk) = e3t_0(jk)
873	e3v (:,:,jk) = e3t_0(jk)
874	e3uw(:,:,jk) = e3w_0(jk)
875	e3vw(:,:,jk) = e3w_0(jk)
876	END DO
877	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
878	DO jj = 1, jpjm1
879	DO ji = 1, fs_jpim1 ! vector opt.
880	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
881	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
882	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
883	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
884	END DO
885	END DO
886	END DO
887	! ! Boundary conditions
888	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
889	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
890	!
891	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
892	WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk)
893	WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk)
894	WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk)
895	WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk)
896	END DO
897
898	! Scale factor at F-point
899	DO jk = 1, jpk ! initialisation to z-scale factors
900	e3f (:,:,jk) = e3t_0(jk)
901	END DO
902	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
903	DO jj = 1, jpjm1
904	DO ji = 1, fs_jpim1 ! vector opt.
905	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
906	END DO
907	END DO
908	END DO
909	CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions
910	!
911	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
912	WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk)
913	END DO
914	!!gm bug ? : must be a do loop with mj0,mj1
915	!
916	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
917	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
918	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
919	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
920	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
921
922	! Control of the sign
923	IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
924	IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
925	IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
926	IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
927
928	! Compute gdep3w (vertical sum of e3w)
929	gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
930	DO jk = 2, jpk
931	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
932	END DO
933
934	! ! ================= !
935	IF(lwp .AND. ll_print) THEN ! Control print !
936	! ! ================= !
937	DO jj = 1,jpj
938	DO ji = 1, jpi
939	ik = MAX( mbathy(ji,jj), 1 )
940	zprt(ji,jj,1) = e3t (ji,jj,ik)
941	zprt(ji,jj,2) = e3w (ji,jj,ik)
942	zprt(ji,jj,3) = e3u (ji,jj,ik)
943	zprt(ji,jj,4) = e3v (ji,jj,ik)
944	zprt(ji,jj,5) = e3f (ji,jj,ik)
945	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
946	END DO
947	END DO
948	WRITE(numout,*)
949	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
950	WRITE(numout,*)
951	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
952	WRITE(numout,*)
953	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
954	WRITE(numout,*)
955	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
956	WRITE(numout,*)
957	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
958	WRITE(numout,*)
959	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
960	ENDIF
961
962	! ! =============== !
963	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
964	! ! =============== !
965
966	! ! =================== !
967	IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check !
968	! ! =================== !
969	END SUBROUTINE zgr_zps
970
971
972	FUNCTION fssig( pk ) RESULT( pf )
973	!!----------------------------------------------------------------------
974	!! * ROUTINE eos_init *
975	!!
976	!! ** Purpose : provide the analytical function in s-coordinate
977	!!
978	!! ** Method : the function provide the non-dimensional position of
979	!! T and W (i.e. between 0 and 1)
980	!! T-points at integer values (between 1 and jpk)
981	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
982	!!
983	!! Reference : ???
984	!!----------------------------------------------------------------------
985	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
986	REAL(wp) :: pf ! sigma value
987	!!----------------------------------------------------------------------
988	!
989	pf = ( TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb ) ) &
990	& - TANH( thetb * theta ) ) &
991	& * ( COSH( theta ) &
992	& + COSH( theta * ( 2.e0 * thetb - 1.e0 ) ) ) &
993	& / ( 2.e0 * SINH( theta ) )
994	!
995	END FUNCTION fssig
996
997
998	SUBROUTINE zgr_sco
999	!!----------------------------------------------------------------------
1000	!! * ROUTINE zgr_sco *
1001	!!
1002	!! ** Purpose : define the s-coordinate system
1003	!!
1004	!! ** Method : s-coordinate
1005	!! The depth of model levels is defined as the product of an
1006	!! analytical function by the local bathymetry, while the vertical
1007	!! scale factors are defined as the product of the first derivative
1008	!! of the analytical function by the bathymetry.
1009	!! (this solution save memory as depth and scale factors are not
1010	!! 3d fields)
1011	!! - Read bathymetry (in meters) at t-point and compute the
1012	!! bathymetry at u-, v-, and f-points.
1013	!! hbatu = mi( hbatt )
1014	!! hbatv = mj( hbatt )
1015	!! hbatf = mi( mj( hbatt ) )
1016	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
1017	!! function and its derivative given as function.
1018	!! gsigt(k) = fssig (k )
1019	!! gsigw(k) = fssig (k-0.5)
1020	!! esigt(k) = fsdsig(k )
1021	!! esigw(k) = fsdsig(k-0.5)
1022	!! This routine is given as an example, it must be modified
1023	!! following the user s desiderata. nevertheless, the output as
1024	!! well as the way to compute the model levels and scale factors
1025	!! must be respected in order to insure second order a!!uracy
1026	!! schemes.
1027	!!
1028	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1029	!!----------------------------------------------------------------------
1030	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1031	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1032	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars
1033	REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace
1034	REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - -
1035	!!
1036	NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
1037	!!----------------------------------------------------------------------
1038
1039	REWIND( numnam ) ! Read Namelist nam_zgr_sco : sigma-stretching parameters
1040	READ ( numnam, nam_zgr_sco )
1041
1042	IF(lwp) THEN ! control print
1043	WRITE(numout,*)
1044	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1045	WRITE(numout,*) '~~~~~~~~~~~'
1046	WRITE(numout,*) ' Namelist nam_zgr_sco'
1047	WRITE(numout,*) ' sigma-stretching coeffs '
1048	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max
1049	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min
1050	WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta
1051	WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb
1052	WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max
1053	ENDIF
1054
1055	hift(:,:) = sbot_min ! set the minimum depth for the s-coordinate
1056	hifu(:,:) = sbot_min
1057	hifv(:,:) = sbot_min
1058	hiff(:,:) = sbot_min
1059	! ! set maximum ocean depth
1060	bathy(:,:) = MIN( sbot_max, bathy(:,:) )
1061
1062	! ! =============================
1063	! ! Define the envelop bathymetry (hbatt)
1064	! ! =============================
1065	! Smooth the bathymetry (if required)
1066	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
1067	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1068	!
1069	! use r-value to create hybrid coordinates
1070	DO jj = 1, jpj
1071	DO ji = 1, jpi
1072	zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
1073	END DO
1074	END DO
1075	jl = 0
1076	zrmax = 1.e0
1077	! ! ================ !
1078	DO WHILE ( jl <= 10000 .AND. zrmax > r_max ) ! Iterative loop !
1079	! ! ================ !
1080	jl = jl + 1
1081	zrmax = 0.e0
1082	zmsk(:,:) = 0.e0
1083	DO jj = 1, nlcj
1084	DO ji = 1, nlci
1085	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1086	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1087	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1088	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1089	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1090	IF( zri(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1091	IF( zri(ji,jj) > r_max ) zmsk(iip1,jj ) = 1.0
1092	IF( zrj(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1093	IF( zrj(ji,jj) > r_max ) zmsk(ji ,ijp1) = 1.0
1094	END DO
1095	END DO
1096	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1097	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1098	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
1099	DO jj = 1, nlcj
1100	DO ji = 1, nlci
1101	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1102	END DO
1103	END DO
1104	!
1105	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1106	!
1107	DO jj = 1, nlcj
1108	DO ji = 1, nlci
1109	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1110	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1111	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1112	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1113	IF( zmsk(ji,jj) == 1.0 ) THEN
1114	ztmp(ji,jj) = ( &
1115	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1116	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
1117	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1118	& ) / ( &
1119	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1120	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
1121	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1122	& )
1123	ENDIF
1124	END DO
1125	END DO
1126	!
1127	DO jj = 1, nlcj
1128	DO ji = 1, nlci
1129	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1130	END DO
1131	END DO
1132	!
1133	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
1134	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
1135	DO jj = 1, nlcj
1136	DO ji = 1, nlci
1137	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
1138	END DO
1139	END DO
1140	! ! ================ !
1141	END DO ! End loop !
1142	! ! ================ !
1143	!
1144	! ! envelop bathymetry saved in hbatt
1145	hbatt(:,:) = zenv(:,:)
1146	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
1147	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1148	DO jj = 1, jpj
1149	DO ji = 1, jpi
1150	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1151	hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1152	END DO
1153	END DO
1154	ENDIF
1155	!
1156	IF(lwp) THEN ! Control print
1157	WRITE(numout,*)
1158	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1159	WRITE(numout,*)
1160	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1161	IF( nprint == 1 ) THEN
1162	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1163	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1164	ENDIF
1165	ENDIF
1166
1167	! ! ==============================
1168	! ! hbatu, hbatv, hbatf fields
1169	! ! ==============================
1170	IF(lwp) THEN
1171	WRITE(numout,*)
1172	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
1173	ENDIF
1174	hbatu(:,:) = sbot_min
1175	hbatv(:,:) = sbot_min
1176	hbatf(:,:) = sbot_min
1177	DO jj = 1, jpjm1
1178	DO ji = 1, jpim1
1179	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1180	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1181	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1182	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1183	END DO
1184	END DO
1185	!
1186	! Apply lateral boundary condition
1187	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1188	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
1189	DO jj = 1, jpj
1190	DO ji = 1, jpi
1191	IF( hbatu(ji,jj) == 0.e0 ) THEN
1192	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min
1193	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1194	ENDIF
1195	END DO
1196	END DO
1197	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
1198	DO jj = 1, jpj
1199	DO ji = 1, jpi
1200	IF( hbatv(ji,jj) == 0.e0 ) THEN
1201	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min
1202	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1203	ENDIF
1204	END DO
1205	END DO
1206	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
1207	DO jj = 1, jpj
1208	DO ji = 1, jpi
1209	IF( hbatf(ji,jj) == 0.e0 ) THEN
1210	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min
1211	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1212	ENDIF
1213	END DO
1214	END DO
1215
1216	!!bug: key_helsinki a verifer
1217	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1218	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1219	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1220	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1221
1222	IF( nprint == 1 .AND. lwp ) THEN
1223	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1224	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1225	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1226	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1227	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1228	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1229	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1230	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1231	ENDIF
1232	!! helsinki
1233
1234	! ! =======================
1235	! ! s-ccordinate fields (gdep., e3.)
1236	! ! =======================
1237	!
1238	! non-dimensional "sigma" for model level depth at w- and t-levels
1239	DO jk = 1, jpk
1240	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1241	gsigt(jk) = -fssig( REAL(jk,wp) )
1242	END DO
1243	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1244	!
1245	! Coefficients for vertical scale factors at w-, t- levels
1246	!!gm bug : define it from analytical function, not like juste bellow....
1247	!!gm or betteroffer the 2 possibilities....
1248	DO jk = 1, jpkm1
1249	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
1250	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1251	END DO
1252	esigw( 1 ) = esigw( 2 )
1253	esigt(jpk) = esigt(jpkm1)
1254
1255	!!gm original form
1256	!!org DO jk = 1, jpk
1257	!!org esigt(jk)=fsdsig( FLOAT(jk) )
1258	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1259	!!org END DO
1260	!!gm
1261	!
1262	! Coefficients for vertical depth as the sum of e3w scale factors
1263	gsi3w(1) = 0.5 * esigw(1)
1264	DO jk = 2, jpk
1265	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1266	END DO
1267	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1268	DO jk = 1, jpk
1269	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1270	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1271	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1272	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1273	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoefw)
1274	END DO
1275	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1276	DO jj = 1, jpj
1277	DO ji = 1, jpi
1278	DO jk = 1, jpk
1279	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1280	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1281	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1282	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1283	!
1284	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1285	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1286	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1287	END DO
1288	END DO
1289	END DO
1290	!
1291	! HYBRID :
1292	DO jj = 1, jpj
1293	DO ji = 1, jpi
1294	DO jk = 1, jpkm1
1295	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1296	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1297	END DO
1298	END DO
1299	END DO
1300	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1301	& ' MAX ', MAXVAL( mbathy(:,:) )
1302
1303
1304	! ! ===========
1305	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
1306	! ! ===========
1307
1308	! ! =============
1309	IF(lwp) THEN ! Control print
1310	! ! =============
1311	WRITE(numout,*)
1312	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1313	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
1314	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1315	ENDIF
1316	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1317	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1318	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
1319	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
1320	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
1321	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
1322	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
1323	& ' w ', MINVAL( fse3w (:,:,:) )
1324
1325	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
1326	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
1327	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
1328	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
1329	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
1330	& ' w ', MAXVAL( fse3w (:,:,:) )
1331	ENDIF
1332	!
1333	IF(lwp) THEN ! selected vertical profiles
1334	WRITE(numout,*)
1335	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1336	WRITE(numout,*) ' ~~~~~~ --------------------'
1337	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1338	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
1339	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1340	DO jj = mj0(20), mj1(20)
1341	DO ji = mi0(20), mi1(20)
1342	WRITE(numout,*)
1343	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1344	WRITE(numout,*) ' ~~~~~~ --------------------'
1345	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1346	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1347	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1348	END DO
1349	END DO
1350	DO jj = mj0(74), mj1(74)
1351	DO ji = mi0(100), mi1(100)
1352	WRITE(numout,*)
1353	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1354	WRITE(numout,*) ' ~~~~~~ --------------------'
1355	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1356	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1357	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1358	END DO
1359	END DO
1360	ENDIF
1361
1362	!!gm bug? no more necessary? if ! defined key_helsinki
1363	DO jk = 1, jpk
1364	DO jj = 1, jpj
1365	DO ji = 1, jpi
1366	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1367	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1368	CALL ctl_stop( ctmp1 )
1369	ENDIF
1370	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1371	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1372	CALL ctl_stop( ctmp1 )
1373	ENDIF
1374	END DO
1375	END DO
1376	END DO
1377	!!gm bug #endif
1378	!
1379	END SUBROUTINE zgr_sco
1380
1381	#endif
1382
1383	!!======================================================================
1384	END MODULE domzgr

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

Download in other formats: