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

Last change on this file since 1099 was 1099, checked in by ctlod, 16 years ago
trunk: code review stuff, see tickets: #205, #191, #130
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 65.3 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
7	!! ! 1997-07 (G. Madec) lbc_lnk call
8	!! ! 1997-04 (J.-O. Beismann)
9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
15	!!----------------------------------------------------------------------
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 :: ii_bump, ij_bump, ih ! bump center position
267	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
268	REAL(wp) :: zi , zj , zh ! temporary scalars
269	INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data
270	REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data
271	!!----------------------------------------------------------------------
272
273	IF(lwp) WRITE(numout,*)
274	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
275	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
276
277	! ! ================== !
278	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
279	! ! ================== !
280	! ! global domain level and meter bathymetry (idta,zdta)
281	!
282	IF( ntopo == 0 ) THEN ! flat basin
283	IF(lwp) WRITE(numout,*)
284	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
285	idta(:,:) = jpkm1 ! before last level
286	zdta(:,:) = gdepw_0(jpk) ! last w-point depth
287	h_oce = gdepw_0(jpk)
288	ELSE ! bump centered in the basin
289	IF(lwp) WRITE(numout,*)
290	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
291	ii_bump = jpidta / 2 ! i-index of the bump center
292	ij_bump = jpjdta / 2 ! j-index of the bump center
293	r_bump = 50000.e0 ! bump radius (meters)
294	h_bump = 2700.e0 ! bump height (meters)
295	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
296	IF(lwp) WRITE(numout,*) ' bump characteristics: '
297	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
298	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
299	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
300	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
301	!
302	DO jj = 1, jpjdta ! zdta :
303	DO ji = 1, jpidta
304	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
305	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
306	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
307	END DO
308	END DO
309	! ! idta :
310	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
311	idta(:,:) = jpkm1
312	ELSE ! z-coordinate (zco or zps): step-like topography
313	idta(:,:) = jpkm1
314	DO jk = 1, jpkm1
315	DO jj = 1, jpjdta
316	DO ji = 1, jpidta
317	IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk
318	END DO
319	END DO
320	END DO
321	ENDIF
322	ENDIF
323	! ! set GLOBAL boundary conditions
324	! ! Caution : idta on the global domain: use of jperio, not nperio
325	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
326	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0
327	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
328	ELSEIF( jperio == 2 ) THEN
329	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
330	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
331	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0
332	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0
333	ELSE
334	ih = 0 ; zh = 0.e0
335	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
336	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
337	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
338	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
339	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
340	ENDIF
341
342	! ! local domain level and meter bathymetries (mbathy,bathy)
343	mbathy(:,:) = 0 ! set to zero extra halo points
344	bathy (:,:) = 0.e0 ! (require for mpp case)
345	DO jj = 1, nlcj ! interior values
346	DO ji = 1, nlci
347	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
348	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
349	END DO
350	END DO
351	!
352	! ! ================ !
353	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
354	! ! ================ !
355	!
356	IF( ln_zco ) THEN ! zco : read level bathymetry
357	CALL iom_open( 'bathy_level.nc', inum )
358	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
359	CALL iom_close (inum)
360	mbathy(:,:) = INT( bathy(:,:) )
361	ENDIF
362	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
363	CALL iom_open( 'bathy_meter.nc', inum )
364	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
365	CALL iom_close (inum)
366	ENDIF
367	! ! =============== !
368	ELSE ! error !
369	! ! =============== !
370	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
371	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
372	ENDIF
373	!
374	! ! =========================== !
375	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
376	DO jl = 1, jpncs ! =========================== !
377	DO jj = ncsj1(jl), ncsj2(jl)
378	DO ji = ncsi1(jl), ncsi2(jl)
379	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
380	bathy (ji,jj) = 0.e0
381	END DO
382	END DO
383	END DO
384	ENDIF
385	#if defined key_orca_lev10
386	! ! ================= !
387	mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !
388	! ! ================= !
389	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
390	#endif
391	! ! =============== !
392	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
393	! ! =============== !
394
395	! ! =================== !
396	IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check !
397	! ! =================== !
398	END SUBROUTINE zgr_bat
399
400
401	SUBROUTINE zgr_bat_zoom
402	!!----------------------------------------------------------------------
403	!! * ROUTINE zgr_bat_zoom *
404	!!
405	!! ** Purpose : - Close zoom domain boundary if necessary
406	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
407	!!
408	!! ** Method :
409	!!
410	!! ** Action : - update mbathy: level bathymetry (in level index)
411	!!----------------------------------------------------------------------
412	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
413	!!----------------------------------------------------------------------
414	!
415	IF(lwp) WRITE(numout,*)
416	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
417	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
418	!
419	! Zoom domain
420	! ===========
421	!
422	! Forced closed boundary if required
423	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
424	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
425	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
426	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
427	!
428	! Configuration specific domain modifications
429	! (here, ORCA arctic configuration: suppress Med Sea)
430	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
431	SELECT CASE ( jp_cfg )
432	! ! =======================
433	CASE ( 2 ) ! ORCA_R2 configuration
434	! ! =======================
435	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
436	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
437	ij0 = 98 ; ij1 = 110
438	! ! =======================
439	CASE ( 05 ) ! ORCA_R05 configuration
440	! ! =======================
441	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
442	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
443	ij0 = 314 ; ij1 = 370
444	END SELECT
445	!
446	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
447	!
448	ENDIF
449	!
450	END SUBROUTINE zgr_bat_zoom
451
452
453	SUBROUTINE zgr_bat_ctl
454	!!----------------------------------------------------------------------
455	!! * ROUTINE zgr_bat_ctl *
456	!!
457	!! ** Purpose : check the bathymetry in levels
458	!!
459	!! ** Method : The array mbathy is checked to verified its consistency
460	!! with the model options. in particular:
461	!! mbathy must have at least 1 land grid-points (mbathy<=0)
462	!! along closed boundary.
463	!! mbathy must be cyclic IF jperio=1.
464	!! mbathy must be lower or equal to jpk-1.
465	!! isolated ocean grid points are suppressed from mbathy
466	!! since they are only connected to remaining
467	!! ocean through vertical diffusion.
468	!! C A U T I O N : mbathy will be modified during the initializa-
469	!! tion phase to become the number of non-zero w-levels of a water
470	!! column, with a minimum value of 1.
471	!!
472	!! ** Action : - update mbathy: level bathymetry (in level index)
473	!! - update bathy : meter bathymetry (in meters)
474	!!----------------------------------------------------------------------
475	INTEGER :: ji, jj, jl ! dummy loop indices
476	INTEGER :: icompt, ibtest, ikmax ! temporary integers
477	REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace
478	!!----------------------------------------------------------------------
479
480	IF(lwp) WRITE(numout,*)
481	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
482	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
483
484	! ! Suppress isolated ocean grid points
485	IF(lwp) WRITE(numout,*)
486	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
487	IF(lwp) WRITE(numout,*)' -----------------------------------'
488	icompt = 0
489	DO jl = 1, 2
490	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
491	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
492	mbathy(jpi,:) = mbathy( 2 ,:)
493	ENDIF
494	DO jj = 2, jpjm1
495	DO ji = 2, jpim1
496	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
497	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
498	IF( ibtest < mbathy(ji,jj) ) THEN
499	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
500	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
501	mbathy(ji,jj) = ibtest
502	icompt = icompt + 1
503	ENDIF
504	END DO
505	END DO
506	END DO
507	IF( icompt == 0 ) THEN
508	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
509	ELSE
510	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
511	ENDIF
512	IF( lk_mpp ) THEN
513	zbathy(:,:) = FLOAT( mbathy(:,:) )
514	CALL lbc_lnk( zbathy, 'T', 1. )
515	mbathy(:,:) = INT( zbathy(:,:) )
516	ENDIF
517
518	! ! East-west cyclic boundary conditions
519	IF( nperio == 0 ) THEN
520	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
521	IF( lk_mpp ) THEN
522	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
523	IF( jperio /= 1 ) mbathy(1,:) = 0
524	ENDIF
525	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
526	IF( jperio /= 1 ) mbathy(nlci,:) = 0
527	ENDIF
528	ELSE
529	IF( ln_zco .OR. ln_zps ) THEN
530	mbathy( 1 ,:) = 0
531	mbathy(jpi,:) = 0
532	ELSE
533	mbathy( 1 ,:) = jpkm1
534	mbathy(jpi,:) = jpkm1
535	ENDIF
536	ENDIF
537	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
538	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
539	mbathy( 1 ,:) = mbathy(jpim1,:)
540	mbathy(jpi,:) = mbathy( 2 ,:)
541	ELSEIF( nperio == 2 ) THEN
542	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
543	ELSE
544	IF(lwp) WRITE(numout,*) ' e r r o r'
545	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
546	! STOP 'dom_mba'
547	ENDIF
548
549	! Set to zero mbathy over islands if necessary (lk_isl=F)
550	IF( .NOT. lk_isl ) THEN ! No island
551	IF(lwp) WRITE(numout,*)
552	IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'
553	IF(lwp) WRITE(numout,*) ' ----------------------------'
554	!
555	mbathy(:,:) = MAX( 0, mbathy(:,:) )
556	!
557	! Boundary condition on mbathy
558	IF( .NOT.lk_mpp ) THEN
559	!!gm !!bug ??? think about it !
560	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
561	zbathy(:,:) = FLOAT( mbathy(:,:) )
562	CALL lbc_lnk( zbathy, 'T', 1. )
563	mbathy(:,:) = INT( zbathy(:,:) )
564	ENDIF
565	ENDIF
566
567	! Number of ocean level inferior or equal to jpkm1
568	ikmax = 0
569	DO jj = 1, jpj
570	DO ji = 1, jpi
571	ikmax = MAX( ikmax, mbathy(ji,jj) )
572	END DO
573	END DO
574	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
575	IF( ikmax > jpkm1 ) THEN
576	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
577	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
578	ELSE IF( ikmax < jpkm1 ) THEN
579	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
580	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
581	ENDIF
582
583	! control print
584	IF( lwp .AND. nprint == 1 ) THEN
585	WRITE(numout,*)
586	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
587	WRITE(numout,*) ' ------------------'
588	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
589	WRITE(numout,*)
590	ENDIF
591	!
592	END SUBROUTINE zgr_bat_ctl
593
594
595	SUBROUTINE zgr_zco
596	!!----------------------------------------------------------------------
597	!! * ROUTINE zgr_zco *
598	!!
599	!! ** Purpose : define the z-coordinate system
600	!!
601	!! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T
602	!!----------------------------------------------------------------------
603	INTEGER :: jk
604	!!----------------------------------------------------------------------
605	!
606	IF( .NOT.lk_zco ) THEN
607	DO jk = 1, jpk
608	fsdept(:,:,jk) = gdept_0(jk)
609	fsdepw(:,:,jk) = gdepw_0(jk)
610	fsde3w(:,:,jk) = gdepw_0(jk)
611	fse3t (:,:,jk) = e3t_0(jk)
612	fse3u (:,:,jk) = e3t_0(jk)
613	fse3v (:,:,jk) = e3t_0(jk)
614	fse3f (:,:,jk) = e3t_0(jk)
615	fse3w (:,:,jk) = e3w_0(jk)
616	fse3uw(:,:,jk) = e3w_0(jk)
617	fse3vw(:,:,jk) = e3w_0(jk)
618	END DO
619	ENDIF
620	!
621	END SUBROUTINE zgr_zco
622
623	#if defined key_zco
624	!!----------------------------------------------------------------------
625	!! 'key_zco' : "pure" zco (gdep & e3 are 1D arrays)
626	!!----------------------------------------------------------------------
627	SUBROUTINE zgr_zps ! Empty routine
628	END SUBROUTINE zgr_zps
629	SUBROUTINE zgr_sco ! Empty routine
630	END SUBROUTINE zgr_sco
631
632	#else
633	!!----------------------------------------------------------------------
634	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
635	!!----------------------------------------------------------------------
636
637	SUBROUTINE zgr_zps
638	!!----------------------------------------------------------------------
639	!! * ROUTINE zgr_zps *
640	!!
641	!! ** Purpose : the depth and vertical scale factor in partial step
642	!! z-coordinate case
643	!!
644	!! ** Method : Partial steps : computes the 3D vertical scale factors
645	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
646	!! a partial step representation of bottom topography.
647	!!
648	!! The reference depth of model levels is defined from an analytical
649	!! function the derivative of which gives the reference vertical
650	!! scale factors.
651	!! From depth and scale factors reference, we compute there new value
652	!! with partial steps on 3d arrays ( i, j, k ).
653	!!
654	!! w-level: gdepw(i,j,k) = fsdep(k)
655	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
656	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
657	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
658	!!
659	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
660	!! we find the mbathy index of the depth at each grid point.
661	!! This leads us to three cases:
662	!!
663	!! - bathy = 0 => mbathy = 0
664	!! - 1 < mbathy < jpkm1
665	!! - bathy > gdepw(jpk) => mbathy = jpkm1
666	!!
667	!! Then, for each case, we find the new depth at t- and w- levels
668	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
669	!! and f-points.
670	!!
671	!! This routine is given as an example, it must be modified
672	!! following the user s desiderata. nevertheless, the output as
673	!! well as the way to compute the model levels and scale factors
674	!! must be respected in order to insure second order accuracy
675	!! schemes.
676	!!
677	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
678	!! - - - - - - - gdept, gdepw and e3. are positives
679	!!
680	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
681	!!----------------------------------------------------------------------
682	INTEGER :: ji, jj, jk ! dummy loop indices
683	INTEGER :: ik, it ! temporary integers
684	LOGICAL :: ll_print ! Allow control print for debugging
685	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
686	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
687	REAL(wp) :: zmax, zmin ! Maximum and minimum depth
688	REAL(wp) :: zdiff ! temporary scalar
689	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace
690	!!---------------------------------------------------------------------
691
692	IF(lwp) WRITE(numout,*)
693	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
694	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
695	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
696
697	ll_print=.FALSE. ! Local variable for debugging
698	!! ll_print=.TRUE.
699
700	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
701	WRITE(numout,*)
702	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
703	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
704	ENDIF
705
706
707	! bathymetry in level (from bathy_meter)
708	! ===================
709	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
710	zmin = gdepw_0(4) ! minimum depth = 3 levels
711
712	mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available
713
714	! ! storage of land and island's number (zera and negative values) in mbathy
715	WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )
716
717	! bounded value of bathy
718	!!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0
719	DO jj = 1, jpj
720	DO ji= 1, jpi
721	IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0
722	ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )
723	ENDIF
724	END DO
725	END DO
726
727	! Compute mbathy for ocean points (i.e. the number of ocean levels)
728	! find the number of ocean levels such that the last level thickness
729	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
730	! e3t_0 is the reference level thickness
731	DO jk = jpkm1, 1, -1
732	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
733	DO jj = 1, jpj
734	DO ji = 1, jpi
735	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
736	END DO
737	END DO
738	END DO
739
740	! Scale factors and depth at T- and W-points
741	DO jk = 1, jpk ! intitialization to the reference z-coordinate
742	gdept(:,:,jk) = gdept_0(jk)
743	gdepw(:,:,jk) = gdepw_0(jk)
744	e3t (:,:,jk) = e3t_0 (jk)
745	e3w (:,:,jk) = e3w_0 (jk)
746	END DO
747	hdept(:,:) = gdept(:,:,2 )
748	hdepw(:,:) = gdepw(:,:,3 )
749	!
750	DO jj = 1, jpj
751	DO ji = 1, jpi
752	ik = mbathy(ji,jj)
753	IF( ik > 0 ) THEN ! ocean point only
754	! max ocean level case
755	IF( ik == jpkm1 ) THEN
756	zdepwp = bathy(ji,jj)
757	ze3tp = bathy(ji,jj) - gdepw_0(ik)
758	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
759	e3t(ji,jj,ik ) = ze3tp
760	e3t(ji,jj,ik+1) = ze3tp
761	e3w(ji,jj,ik ) = ze3wp
762	e3w(ji,jj,ik+1) = ze3tp
763	gdepw(ji,jj,ik+1) = zdepwp
764	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
765	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
766	!
767	ELSE ! standard case
768	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
769	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
770	ENDIF
771	!gm Bug? check the gdepw_0
772	! ... on ik
773	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
774	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
775	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
776	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
777	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
778	e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
779	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
780	! ... on ik+1
781	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
782	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
783	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
784	ENDIF
785	ENDIF
786	END DO
787	END DO
788	!
789	it = 0
790	DO jj = 1, jpj
791	DO ji = 1, jpi
792	ik = mbathy(ji,jj)
793	IF( ik > 0 ) THEN ! ocean point only
794	hdept(ji,jj) = gdept(ji,jj,ik )
795	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
796	e3tp (ji,jj) = e3t(ji,jj,ik )
797	e3wp (ji,jj) = e3w(ji,jj,ik )
798	! test
799	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
800	IF( zdiff <= 0. .AND. lwp ) THEN
801	it = it + 1
802	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
803	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
804	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
805	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
806	ENDIF
807	ENDIF
808	END DO
809	END DO
810
811	! Scale factors and depth at U-, V-, UW and VW-points
812	DO jk = 1, jpk ! initialisation to z-scale factors
813	e3u (:,:,jk) = e3t_0(jk)
814	e3v (:,:,jk) = e3t_0(jk)
815	e3uw(:,:,jk) = e3w_0(jk)
816	e3vw(:,:,jk) = e3w_0(jk)
817	END DO
818	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
819	DO jj = 1, jpjm1
820	DO ji = 1, fs_jpim1 ! vector opt.
821	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
822	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
823	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
824	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
825	END DO
826	END DO
827	END DO
828	! ! Boundary conditions
829	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
830	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
831	!
832	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
833	WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk)
834	WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk)
835	WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk)
836	WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk)
837	END DO
838
839	! Scale factor at F-point
840	DO jk = 1, jpk ! initialisation to z-scale factors
841	e3f (:,:,jk) = e3t_0(jk)
842	END DO
843	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
844	DO jj = 1, jpjm1
845	DO ji = 1, fs_jpim1 ! vector opt.
846	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
847	END DO
848	END DO
849	END DO
850	CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions
851	!
852	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
853	WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk)
854	END DO
855	!!gm bug ? : must be a do loop with mj0,mj1
856	!
857	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
858	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
859	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
860	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
861	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
862
863	! Control of the sign
864	IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
865	IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
866	IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
867	IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
868
869	! Compute gdep3w (vertical sum of e3w)
870	gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
871	DO jk = 2, jpk
872	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
873	END DO
874
875	! ! ================= !
876	IF(lwp .AND. ll_print) THEN ! Control print !
877	! ! ================= !
878	DO jj = 1,jpj
879	DO ji = 1, jpi
880	ik = MAX( mbathy(ji,jj), 1 )
881	zprt(ji,jj,1) = e3t (ji,jj,ik)
882	zprt(ji,jj,2) = e3w (ji,jj,ik)
883	zprt(ji,jj,3) = e3u (ji,jj,ik)
884	zprt(ji,jj,4) = e3v (ji,jj,ik)
885	zprt(ji,jj,5) = e3f (ji,jj,ik)
886	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
887	END DO
888	END DO
889	WRITE(numout,*)
890	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
891	WRITE(numout,*)
892	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
893	WRITE(numout,*)
894	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
895	WRITE(numout,*)
896	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
897	WRITE(numout,*)
898	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
899	WRITE(numout,*)
900	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
901	ENDIF
902
903	! ! =============== !
904	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
905	! ! =============== !
906
907	! ! =================== !
908	IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check !
909	! ! =================== !
910	END SUBROUTINE zgr_zps
911
912
913	FUNCTION fssig( pk ) RESULT( pf )
914	!!----------------------------------------------------------------------
915	!! * ROUTINE eos_init *
916	!!
917	!! ** Purpose : provide the analytical function in s-coordinate
918	!!
919	!! ** Method : the function provide the non-dimensional position of
920	!! T and W (i.e. between 0 and 1)
921	!! T-points at integer values (between 1 and jpk)
922	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
923	!!
924	!! Reference : ???
925	!!----------------------------------------------------------------------
926	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
927	REAL(wp) :: pf ! sigma value
928	!!----------------------------------------------------------------------
929	!
930	pf = ( TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb ) ) &
931	& - TANH( thetb * theta ) ) &
932	& * ( COSH( theta ) &
933	& + COSH( theta * ( 2.e0 * thetb - 1.e0 ) ) ) &
934	& / ( 2.e0 * SINH( theta ) )
935	!
936	END FUNCTION fssig
937
938
939	SUBROUTINE zgr_sco
940	!!----------------------------------------------------------------------
941	!! * ROUTINE zgr_sco *
942	!!
943	!! ** Purpose : define the s-coordinate system
944	!!
945	!! ** Method : s-coordinate
946	!! The depth of model levels is defined as the product of an
947	!! analytical function by the local bathymetry, while the vertical
948	!! scale factors are defined as the product of the first derivative
949	!! of the analytical function by the bathymetry.
950	!! (this solution save memory as depth and scale factors are not
951	!! 3d fields)
952	!! - Read bathymetry (in meters) at t-point and compute the
953	!! bathymetry at u-, v-, and f-points.
954	!! hbatu = mi( hbatt )
955	!! hbatv = mj( hbatt )
956	!! hbatf = mi( mj( hbatt ) )
957	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
958	!! function and its derivative given as function.
959	!! gsigt(k) = fssig (k )
960	!! gsigw(k) = fssig (k-0.5)
961	!! esigt(k) = fsdsig(k )
962	!! esigw(k) = fsdsig(k-0.5)
963	!! This routine is given as an example, it must be modified
964	!! following the user s desiderata. nevertheless, the output as
965	!! well as the way to compute the model levels and scale factors
966	!! must be respected in order to insure second order a!!uracy
967	!! schemes.
968	!!
969	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
970	!!----------------------------------------------------------------------
971	INTEGER :: ji, jj, jk, jl ! dummy loop argument
972	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
973	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars
974	REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace
975	REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - -
976	!!
977	NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
978	!!----------------------------------------------------------------------
979
980	REWIND( numnam ) ! Read Namelist nam_zgr_sco : sigma-stretching parameters
981	READ ( numnam, nam_zgr_sco )
982
983	IF(lwp) THEN ! control print
984	WRITE(numout,*)
985	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
986	WRITE(numout,*) '~~~~~~~~~~~'
987	WRITE(numout,*) ' Namelist nam_zgr_sco'
988	WRITE(numout,*) ' sigma-stretching coeffs '
989	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max
990	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min
991	WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta
992	WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb
993	WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max
994	ENDIF
995
996	hift(:,:) = sbot_min ! set the minimum depth for the s-coordinate
997	hifu(:,:) = sbot_min
998	hifv(:,:) = sbot_min
999	hiff(:,:) = sbot_min
1000	! ! set maximum ocean depth
1001	bathy(:,:) = MIN( sbot_max, bathy(:,:) )
1002
1003	! ! =============================
1004	! ! Define the envelop bathymetry (hbatt)
1005	! ! =============================
1006	! Smooth the bathymetry (if required)
1007	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
1008	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1009	!
1010	! use r-value to create hybrid coordinates
1011	DO jj = 1, jpj
1012	DO ji = 1, jpi
1013	zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
1014	END DO
1015	END DO
1016	jl = 0
1017	zrmax = 1.e0
1018	! ! ================ !
1019	DO WHILE ( jl <= 10000 .AND. zrmax > r_max ) ! Iterative loop !
1020	! ! ================ !
1021	jl = jl + 1
1022	zrmax = 0.e0
1023	zmsk(:,:) = 0.e0
1024	DO jj = 1, nlcj
1025	DO ji = 1, nlci
1026	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1027	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1028	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1029	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1030	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1031	IF( zri(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1032	IF( zri(ji,jj) > r_max ) zmsk(iip1,jj ) = 1.0
1033	IF( zrj(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1034	IF( zrj(ji,jj) > r_max ) zmsk(ji ,ijp1) = 1.0
1035	END DO
1036	END DO
1037	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1038	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1039	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
1040	DO jj = 1, nlcj
1041	DO ji = 1, nlci
1042	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1043	END DO
1044	END DO
1045	!
1046	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1047	!
1048	DO jj = 1, nlcj
1049	DO ji = 1, nlci
1050	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1051	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1052	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1053	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1054	IF( zmsk(ji,jj) == 1.0 ) THEN
1055	ztmp(ji,jj) = ( &
1056	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1057	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
1058	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1059	& ) / ( &
1060	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1061	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
1062	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1063	& )
1064	ENDIF
1065	END DO
1066	END DO
1067	!
1068	DO jj = 1, nlcj
1069	DO ji = 1, nlci
1070	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1071	END DO
1072	END DO
1073	!
1074	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
1075	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
1076	DO jj = 1, nlcj
1077	DO ji = 1, nlci
1078	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
1079	END DO
1080	END DO
1081	! ! ================ !
1082	END DO ! End loop !
1083	! ! ================ !
1084	!
1085	! ! envelop bathymetry saved in hbatt
1086	hbatt(:,:) = zenv(:,:)
1087	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
1088	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1089	DO jj = 1, jpj
1090	DO ji = 1, jpi
1091	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1092	hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1093	END DO
1094	END DO
1095	ENDIF
1096	!
1097	IF(lwp) THEN ! Control print
1098	WRITE(numout,*)
1099	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1100	WRITE(numout,*)
1101	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1102	IF( nprint == 1 ) THEN
1103	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1104	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1105	ENDIF
1106	ENDIF
1107
1108	! ! ==============================
1109	! ! hbatu, hbatv, hbatf fields
1110	! ! ==============================
1111	IF(lwp) THEN
1112	WRITE(numout,*)
1113	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
1114	ENDIF
1115	hbatu(:,:) = sbot_min
1116	hbatv(:,:) = sbot_min
1117	hbatf(:,:) = sbot_min
1118	DO jj = 1, jpjm1
1119	DO ji = 1, jpim1
1120	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1121	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1122	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1123	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1124	END DO
1125	END DO
1126	!
1127	! Apply lateral boundary condition
1128	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1129	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
1130	DO jj = 1, jpj
1131	DO ji = 1, jpi
1132	IF( hbatu(ji,jj) == 0.e0 ) THEN
1133	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min
1134	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1135	ENDIF
1136	END DO
1137	END DO
1138	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
1139	DO jj = 1, jpj
1140	DO ji = 1, jpi
1141	IF( hbatv(ji,jj) == 0.e0 ) THEN
1142	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min
1143	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1144	ENDIF
1145	END DO
1146	END DO
1147	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
1148	DO jj = 1, jpj
1149	DO ji = 1, jpi
1150	IF( hbatf(ji,jj) == 0.e0 ) THEN
1151	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min
1152	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1153	ENDIF
1154	END DO
1155	END DO
1156
1157	!!bug: key_helsinki a verifer
1158	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1159	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1160	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1161	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1162
1163	IF( nprint == 1 .AND. lwp ) THEN
1164	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1165	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1166	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1167	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1168	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1169	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1170	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1171	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1172	ENDIF
1173	!! helsinki
1174
1175	! ! =======================
1176	! ! s-ccordinate fields (gdep., e3.)
1177	! ! =======================
1178	!
1179	! non-dimensional "sigma" for model level depth at w- and t-levels
1180	DO jk = 1, jpk
1181	gsigw(jk) = -fssig( FLOAT(jk)-0.5 )
1182	gsigt(jk) = -fssig( FLOAT(jk) )
1183	END DO
1184	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1185	!
1186	! Coefficients for vertical scale factors at w-, t- levels
1187	!!gm bug : define it from analytical function, not like juste bellow....
1188	!!gm or betteroffer the 2 possibilities....
1189	DO jk = 1, jpkm1
1190	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
1191	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1192	END DO
1193	esigw( 1 ) = esigw( 2 )
1194	esigt(jpk) = esigt(jpkm1)
1195
1196	!!gm original form
1197	!!org DO jk = 1, jpk
1198	!!org esigt(jk)=fsdsig( FLOAT(jk) )
1199	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1200	!!org END DO
1201	!!gm
1202	!
1203	! Coefficients for vertical depth as the sum of e3w scale factors
1204	gsi3w(1) = 0.5 * esigw(1)
1205	DO jk = 2, jpk
1206	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1207	END DO
1208	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1209	DO jk = 1, jpk
1210	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1211	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1212	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1213	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1214	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoefw)
1215	END DO
1216	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1217	DO jj = 1, jpj
1218	DO ji = 1, jpi
1219	DO jk = 1, jpk
1220	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1221	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1222	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1223	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1224	!
1225	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1226	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1227	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1228	END DO
1229	END DO
1230	END DO
1231	!
1232	! HYBRID :
1233	DO jj = 1, jpj
1234	DO ji = 1, jpi
1235	DO jk = 1, jpkm1
1236	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1237	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1238	END DO
1239	END DO
1240	END DO
1241	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1242	& ' MAX ', MAXVAL( mbathy(:,:) )
1243
1244
1245	! ! ===========
1246	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
1247	! ! ===========
1248
1249	! ! =============
1250	IF(lwp) THEN ! Control print
1251	! ! =============
1252	WRITE(numout,*)
1253	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1254	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
1255	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1256	ENDIF
1257	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1258	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1259	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
1260	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
1261	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
1262	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
1263	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
1264	& ' w ', MINVAL( fse3w (:,:,:) )
1265
1266	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
1267	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
1268	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
1269	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
1270	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
1271	& ' w ', MAXVAL( fse3w (:,:,:) )
1272	ENDIF
1273	!
1274	IF(lwp) THEN ! selected vertical profiles
1275	WRITE(numout,*)
1276	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1277	WRITE(numout,*) ' ~~~~~~ --------------------'
1278	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1279	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
1280	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1281	DO jj = mj0(20), mj1(20)
1282	DO ji = mi0(20), mi1(20)
1283	WRITE(numout,*)
1284	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1285	WRITE(numout,*) ' ~~~~~~ --------------------'
1286	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1287	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1288	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1289	END DO
1290	END DO
1291	DO jj = mj0(74), mj1(74)
1292	DO ji = mi0(100), mi1(100)
1293	WRITE(numout,*)
1294	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1295	WRITE(numout,*) ' ~~~~~~ --------------------'
1296	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1297	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1298	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1299	END DO
1300	END DO
1301	ENDIF
1302
1303	!!gm bug? no more necessary? if ! defined key_helsinki
1304	DO jk = 1, jpk
1305	DO jj = 1, jpj
1306	DO ji = 1, jpi
1307	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1308	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1309	CALL ctl_stop( ctmp1 )
1310	ENDIF
1311	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1312	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1313	CALL ctl_stop( ctmp1 )
1314	ENDIF
1315	END DO
1316	END DO
1317	END DO
1318	!!gm bug #endif
1319	!
1320	END SUBROUTINE zgr_sco
1321
1322	#endif
1323
1324	!!======================================================================
1325	END MODULE domzgr

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

Download in other formats: