New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domzgr.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 83.8 KB

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

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

Download in other formats: