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 @ 2683

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

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

Download in other formats: