Context Navigation

Back to Ticket #1127

Ticket #1127: domzgr.F90

File domzgr.F90, 81.5 KB (added by poddo, 11 years ago)

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

Download in other formats:

Original Format

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

Context Navigation

Ticket #1127: domzgr.F90

Download in other formats: