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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago
set proper svn properties to all files...
Property svn:keywords set to `Id`
File size: 76.4 KB

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

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

Download in other formats: