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

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3224

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

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

Download in other formats: