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

Context Navigation

source: tags/nemo_v3_2_beta/NEMO/OPA_SRC/DOM/domzgr.F90 @ 8191

Last change on this file since 8191 was 1601, checked in by ctlod, 15 years ago
Doctor naming of OPA namelist variables , see ticket: #526
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 76.9 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$
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	!!gm never do that !!!! ==> Pb at compilation phase on several computer
1086	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3 = 0.0d0
1087	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigt3 = 0.0d0
1088	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsi3w3 = 0.0d0
1089	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigt3 = 0.0d0
1090	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigw3 = 0.0d0
1091	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtu3 = 0.0d0
1092	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtv3 = 0.0d0
1093	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtf3 = 0.0d0
1094	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwu3 = 0.0d0
1095	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwv3 = 0.0d0
1096	!!gm end
1097	!!
1098	NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1099	!!----------------------------------------------------------------------
1100
1101	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
1102	READ ( numnam, namzgr_sco )
1103
1104	IF(lwp) THEN ! control print
1105	WRITE(numout,*)
1106	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1107	WRITE(numout,*) '~~~~~~~~~~~'
1108	WRITE(numout,*) ' Namelist namzgr_sco'
1109	WRITE(numout,*) ' sigma-stretching coeffs '
1110	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max
1111	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min
1112	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta
1113	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb
1114	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax
1115	WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma
1116	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb
1117	WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc
1118	ENDIF
1119
1120	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1121	hifu(:,:) = rn_sbot_min
1122	hifv(:,:) = rn_sbot_min
1123	hiff(:,:) = rn_sbot_min
1124
1125	! ! set maximum ocean depth
1126	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1127
1128	DO jj = 1, jpj
1129	DO ji = 1, jpi
1130	IF (bathy(ji,jj) .gt. 0.0) THEN
1131	bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1132	ENDIF
1133	END DO
1134	END DO
1135	! ! =============================
1136	! ! Define the envelop bathymetry (hbatt)
1137	! ! =============================
1138	! Smooth the bathymetry (if required)
1139	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
1140	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1141	!
1142	! use r-value to create hybrid coordinates
1143	DO jj = 1, jpj
1144	DO ji = 1, jpi
1145	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1146	END DO
1147	END DO
1148	jl = 0
1149	zrmax = 1.e0
1150	! ! ================ !
1151	DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
1152	! ! ================ !
1153	jl = jl + 1
1154	zrmax = 0.e0
1155	zmsk(:,:) = 0.e0
1156	DO jj = 1, nlcj
1157	DO ji = 1, nlci
1158	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1159	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1160	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1161	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1162	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1163	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
1164	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1.0
1165	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
1166	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1.0
1167	END DO
1168	END DO
1169	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1170	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1171	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
1172	DO jj = 1, nlcj
1173	DO ji = 1, nlci
1174	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1175	END DO
1176	END DO
1177	!
1178	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1179	!
1180	DO jj = 1, nlcj
1181	DO ji = 1, nlci
1182	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1183	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1184	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1185	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1186	IF( zmsk(ji,jj) == 1.0 ) THEN
1187	ztmp(ji,jj) = ( &
1188	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1189	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
1190	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1191	& ) / ( &
1192	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1193	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
1194	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1195	& )
1196	ENDIF
1197	END DO
1198	END DO
1199	!
1200	DO jj = 1, nlcj
1201	DO ji = 1, nlci
1202	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1203	END DO
1204	END DO
1205	!
1206	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
1207	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
1208	DO jj = 1, nlcj
1209	DO ji = 1, nlci
1210	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
1211	END DO
1212	END DO
1213	! ! ================ !
1214	END DO ! End loop !
1215	! ! ================ !
1216	!
1217	! ! envelop bathymetry saved in hbatt
1218	hbatt(:,:) = zenv(:,:)
1219	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
1220	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1221	DO jj = 1, jpj
1222	DO ji = 1, jpi
1223	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1224	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1225	END DO
1226	END DO
1227	ENDIF
1228	!
1229	IF(lwp) THEN ! Control print
1230	WRITE(numout,*)
1231	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1232	WRITE(numout,*)
1233	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1234	IF( nprint == 1 ) THEN
1235	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1236	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1237	ENDIF
1238	ENDIF
1239
1240	! ! ==============================
1241	! ! hbatu, hbatv, hbatf fields
1242	! ! ==============================
1243	IF(lwp) THEN
1244	WRITE(numout,*)
1245	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1246	ENDIF
1247	hbatu(:,:) = rn_sbot_min
1248	hbatv(:,:) = rn_sbot_min
1249	hbatf(:,:) = rn_sbot_min
1250	DO jj = 1, jpjm1
1251	DO ji = 1, jpim1
1252	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1253	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1254	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1255	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1256	END DO
1257	END DO
1258	!
1259	! Apply lateral boundary condition
1260	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1261	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
1262	DO jj = 1, jpj
1263	DO ji = 1, jpi
1264	IF( hbatu(ji,jj) == 0.e0 ) THEN
1265	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = rn_sbot_min
1266	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1267	ENDIF
1268	END DO
1269	END DO
1270	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
1271	DO jj = 1, jpj
1272	DO ji = 1, jpi
1273	IF( hbatv(ji,jj) == 0.e0 ) THEN
1274	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = rn_sbot_min
1275	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1276	ENDIF
1277	END DO
1278	END DO
1279	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
1280	DO jj = 1, jpj
1281	DO ji = 1, jpi
1282	IF( hbatf(ji,jj) == 0.e0 ) THEN
1283	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = rn_sbot_min
1284	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1285	ENDIF
1286	END DO
1287	END DO
1288
1289	!!bug: key_helsinki a verifer
1290	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1291	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1292	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1293	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1294
1295	IF( nprint == 1 .AND. lwp ) THEN
1296	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1297	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1298	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1299	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1300	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1301	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1302	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1303	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1304	ENDIF
1305	!! helsinki
1306
1307	! ! =======================
1308	! ! s-ccordinate fields (gdep., e3.)
1309	! ! =======================
1310	!
1311	! non-dimensional "sigma" for model level depth at w- and t-levels
1312
1313	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
1314	! ! below rn_hc, with uniform sigma in shallower waters
1315	DO ji = 1, jpi
1316	DO jj = 1, jpj
1317
1318	IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
1319	DO jk = 1, jpk
1320	gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1321	gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1322	END DO
1323	ELSE ! shallow water, uniform sigma
1324	DO jk = 1, jpk
1325	gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1326	gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1327	END DO
1328	ENDIF
1329	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1330
1331
1332	DO jk = 1, jpkm1
1333	esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1334	esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1335	END DO
1336	esigw3(ji,jj,1) = esigw3(ji,jj,2)
1337	esigt3(ji,jj,jpk) = esigt3(ji,jj,jpkm1)
1338
1339	! Coefficients for vertical depth as the sum of e3w scale factors
1340	gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
1341	DO jk = 2, jpk
1342	gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1343	END DO
1344
1345	DO jk = 1, jpk
1346	zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
1347	zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
1348	gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigt3(ji,jj,jk)+rn_hczcoeft)
1349	gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigw3(ji,jj,jk)+rn_hczcoefw)
1350	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsi3w3(ji,jj,jk)+rn_hczcoefw)
1351	END DO
1352
1353	ENDDO ! for all jj's
1354	ENDDO ! for all ji's
1355
1356
1357	DO ji=1,jpi
1358	DO jj=1,jpj
1359
1360	DO jk = 1, jpk
1361	esigtu3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) ) / &
1362	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1363	esigtv3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji,jj+1)esigt3(ji,jj+1,jk) ) / &
1364	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1365	esigtf3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) + &
1366	hbatt(ji,jj+1)esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)esigt3(ji+1,jj+1,jk) ) / &
1367	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1368	esigwu3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji+1,jj)esigw3(ji+1,jj,jk) ) / &
1369	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1370	esigwv3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji,jj+1)esigw3(ji,jj+1,jk) ) / &
1371	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1372
1373	e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1374	e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1375	e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1376	e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1377	!
1378	e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1379	e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1380	e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1381	END DO
1382
1383	ENDDO
1384	ENDDO
1385
1386	ELSE ! not ln_s_sigma
1387
1388	DO jk = 1, jpk
1389	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1390	gsigt(jk) = -fssig( REAL(jk,wp) )
1391	END DO
1392	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1393	!
1394	! Coefficients for vertical scale factors at w-, t- levels
1395	!!gm bug : define it from analytical function, not like juste bellow....
1396	!!gm or betteroffer the 2 possibilities....
1397	DO jk = 1, jpkm1
1398	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
1399	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1400	END DO
1401	esigw( 1 ) = esigw( 2 )
1402	esigt(jpk) = esigt(jpkm1)
1403
1404	!!gm original form
1405	!!org DO jk = 1, jpk
1406	!!org esigt(jk)=fsdsig( FLOAT(jk) )
1407	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1408	!!org END DO
1409	!!gm
1410	!
1411	! Coefficients for vertical depth as the sum of e3w scale factors
1412	gsi3w(1) = 0.5 * esigw(1)
1413	DO jk = 2, jpk
1414	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1415	END DO
1416	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1417	DO jk = 1, jpk
1418	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1419	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1420	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1421	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1422	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoefw)
1423	END DO
1424	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1425	DO jj = 1, jpj
1426	DO ji = 1, jpi
1427	DO jk = 1, jpk
1428	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1429	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1430	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1431	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1432	!
1433	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1434	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1435	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1436	END DO
1437	END DO
1438	END DO
1439
1440	ENDIF ! ln_s_sigma
1441
1442
1443	!
1444	!! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
1445
1446	fsdept(:,:,:) = gdept (:,:,:)
1447	fsdepw(:,:,:) = gdepw (:,:,:)
1448	fsde3w(:,:,:) = gdep3w(:,:,:)
1449	fse3t (:,:,:) = e3t (:,:,:)
1450	fse3u (:,:,:) = e3u (:,:,:)
1451	fse3v (:,:,:) = e3v (:,:,:)
1452	fse3f (:,:,:) = e3f (:,:,:)
1453	fse3w (:,:,:) = e3w (:,:,:)
1454	fse3uw(:,:,:) = e3uw (:,:,:)
1455	fse3vw(:,:,:) = e3vw (:,:,:)
1456	!!
1457	! HYBRID :
1458	DO jj = 1, jpj
1459	DO ji = 1, jpi
1460	DO jk = 1, jpkm1
1461	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1462	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1463	END DO
1464	END DO
1465	END DO
1466	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1467	& ' MAX ', MAXVAL( mbathy(:,:) )
1468
1469
1470	! ! ===========
1471	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
1472	! ! ===========
1473
1474	! ! =============
1475	IF(lwp) THEN ! Control print
1476	! ! =============
1477	WRITE(numout,*)
1478	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1479	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
1480	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1481	ENDIF
1482	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1483	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1484	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
1485	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
1486	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
1487	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
1488	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
1489	& ' w ', MINVAL( fse3w (:,:,:) )
1490
1491	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
1492	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
1493	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
1494	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
1495	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
1496	& ' w ', MAXVAL( fse3w (:,:,:) )
1497	ENDIF
1498	!
1499	IF(lwp) THEN ! selected vertical profiles
1500	WRITE(numout,*)
1501	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1502	WRITE(numout,*) ' ~~~~~~ --------------------'
1503	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1504	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
1505	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1506	DO jj = mj0(20), mj1(20)
1507	DO ji = mi0(20), mi1(20)
1508	WRITE(numout,*)
1509	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1510	WRITE(numout,*) ' ~~~~~~ --------------------'
1511	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1512	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1513	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1514	END DO
1515	END DO
1516	DO jj = mj0(74), mj1(74)
1517	DO ji = mi0(100), mi1(100)
1518	WRITE(numout,*)
1519	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1520	WRITE(numout,*) ' ~~~~~~ --------------------'
1521	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1522	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1523	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1524	END DO
1525	END DO
1526	ENDIF
1527
1528	!!gm bug? no more necessary? if ! defined key_helsinki
1529	DO jk = 1, jpk
1530	DO jj = 1, jpj
1531	DO ji = 1, jpi
1532	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1533	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1534	CALL ctl_stop( ctmp1 )
1535	ENDIF
1536	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1537	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1538	CALL ctl_stop( ctmp1 )
1539	ENDIF
1540	END DO
1541	END DO
1542	END DO
1543	!!gm bug #endif
1544	!
1545	END SUBROUTINE zgr_sco
1546
1547	#endif
1548
1549	!!======================================================================
1550	END MODULE domzgr

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

Download in other formats: