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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domzgr.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2636

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

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

Download in other formats: