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

Context Navigation

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2236

Last change on this file since 2236 was 2236, checked in by cetlod, 14 years ago
First guess of NEMO_v3.3
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 77.2 KB

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

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

Download in other formats: