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

Context Navigation

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/DOM/domzgr.F90 @ 1878

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

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

Download in other formats: