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

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

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

Context Navigation

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

Last change on this file since 2392 was 2392, checked in by gm, 13 years ago
v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections
Property svn:keywords set to `Id`
File size: 78.9 KB

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

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

Download in other formats: