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

Last change on this file since 2380 was 2380, checked in by acc, 13 years ago
nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)
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( n_cla == 0 ) THEN
415	!
416	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
417	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
418	DO ji = mi0(ii0), mi1(ii1)
419	DO jj = mj0(ij0), mj1(ij1)
420	mbathy(ji,jj) = 15
421	END DO
422	END DO
423	IF(lwp) WRITE(numout,*)
424	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
425	!
426	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
427	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
428	DO ji = mi0(ii0), mi1(ii1)
429	DO jj = mj0(ij0), mj1(ij1)
430	mbathy(ji,jj) = 12
431	END DO
432	END DO
433	IF(lwp) WRITE(numout,*)
434	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
435	!
436	ENDIF
437	!
438	ENDIF
439	!
440	ENDIF
441	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
442	CALL iom_open( 'bathy_meter.nc', inum )
443	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
444	CALL iom_close (inum)
445	! ! =====================
446	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
447	ii0 = 142 ; ii1 = 142 ! Close Halmera Strait
448	ij0 = 51 ; ij1 = 53 ! =====================
449	DO ji = mi0(ii0), mi1(ii1)
450	DO jj = mj0(ij0), mj1(ij1)
451	bathy(ji,jj) = 0.0
452	END DO
453	END DO
454	IF(lwp) WRITE(numout,*)
455	IF(lwp) WRITE(numout,*) ' orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
456	ENDIF
457	! ! =====================
458	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
459	! ! =====================
460	IF( n_cla == 0 ) THEN
461	!
462	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
463	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
464	DO ji = mi0(ii0), mi1(ii1)
465	DO jj = mj0(ij0), mj1(ij1)
466	bathy(ji,jj) = 284.e0
467	END DO
468	END DO
469	IF(lwp) WRITE(numout,*)
470	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
471	!
472	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
473	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
474	DO ji = mi0(ii0), mi1(ii1)
475	DO jj = mj0(ij0), mj1(ij1)
476	bathy(ji,jj) = 137.e0
477	END DO
478	END DO
479	IF(lwp) WRITE(numout,*)
480	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
481	!
482	ENDIF
483	!
484	ENDIF
485	!
486	ENDIF
487	! ! =============== !
488	ELSE ! error !
489	! ! =============== !
490	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
491	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
492	ENDIF
493	!
494	! ! =========================== !
495	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
496	DO jl = 1, jpncs ! =========================== !
497	DO jj = ncsj1(jl), ncsj2(jl)
498	DO ji = ncsi1(jl), ncsi2(jl)
499	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
500	bathy (ji,jj) = 0.e0
501	END DO
502	END DO
503	END DO
504	ENDIF
505	#if defined key_orca_lev10
506	! ! ================= !
507	mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !
508	! ! ================= !
509	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
510	#endif
511	! ! =============== !
512	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
513	! ! =============== !
514
515	#if ! defined key_c1d
516	! ! =================== !
517	CALL zgr_bat_ctl ! Bathymetry check !
518	! ! =================== !
519	#endif
520	END SUBROUTINE zgr_bat
521
522
523	SUBROUTINE zgr_bat_zoom
524	!!----------------------------------------------------------------------
525	!! * ROUTINE zgr_bat_zoom *
526	!!
527	!! ** Purpose : - Close zoom domain boundary if necessary
528	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
529	!!
530	!! ** Method :
531	!!
532	!! ** Action : - update mbathy: level bathymetry (in level index)
533	!!----------------------------------------------------------------------
534	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
535	!!----------------------------------------------------------------------
536	!
537	IF(lwp) WRITE(numout,*)
538	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
539	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
540	!
541	! Zoom domain
542	! ===========
543	!
544	! Forced closed boundary if required
545	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
546	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
547	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
548	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
549	!
550	! Configuration specific domain modifications
551	! (here, ORCA arctic configuration: suppress Med Sea)
552	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
553	SELECT CASE ( jp_cfg )
554	! ! =======================
555	CASE ( 2 ) ! ORCA_R2 configuration
556	! ! =======================
557	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
558	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
559	ij0 = 98 ; ij1 = 110
560	! ! =======================
561	CASE ( 05 ) ! ORCA_R05 configuration
562	! ! =======================
563	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
564	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
565	ij0 = 314 ; ij1 = 370
566	END SELECT
567	!
568	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
569	!
570	ENDIF
571	!
572	END SUBROUTINE zgr_bat_zoom
573
574
575	SUBROUTINE zgr_bat_ctl
576	!!----------------------------------------------------------------------
577	!! * ROUTINE zgr_bat_ctl *
578	!!
579	!! ** Purpose : check the bathymetry in levels
580	!!
581	!! ** Method : The array mbathy is checked to verified its consistency
582	!! with the model options. in particular:
583	!! mbathy must have at least 1 land grid-points (mbathy<=0)
584	!! along closed boundary.
585	!! mbathy must be cyclic IF jperio=1.
586	!! mbathy must be lower or equal to jpk-1.
587	!! isolated ocean grid points are suppressed from mbathy
588	!! since they are only connected to remaining
589	!! ocean through vertical diffusion.
590	!! C A U T I O N : mbathy will be modified during the initializa-
591	!! tion phase to become the number of non-zero w-levels of a water
592	!! column, with a minimum value of 1.
593	!!
594	!! ** Action : - update mbathy: level bathymetry (in level index)
595	!! - update bathy : meter bathymetry (in meters)
596	!!----------------------------------------------------------------------
597	INTEGER :: ji, jj, jl ! dummy loop indices
598	INTEGER :: icompt, ibtest, ikmax ! temporary integers
599	REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace
600	!!----------------------------------------------------------------------
601
602	IF(lwp) WRITE(numout,*)
603	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
604	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
605
606	! ! Suppress isolated ocean grid points
607	IF(lwp) WRITE(numout,*)
608	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
609	IF(lwp) WRITE(numout,*)' -----------------------------------'
610	icompt = 0
611	DO jl = 1, 2
612	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
613	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
614	mbathy(jpi,:) = mbathy( 2 ,:)
615	ENDIF
616	DO jj = 2, jpjm1
617	DO ji = 2, jpim1
618	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
619	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
620	IF( ibtest < mbathy(ji,jj) ) THEN
621	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
622	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
623	mbathy(ji,jj) = ibtest
624	icompt = icompt + 1
625	ENDIF
626	END DO
627	END DO
628	END DO
629	IF( icompt == 0 ) THEN
630	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
631	ELSE
632	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
633	ENDIF
634	IF( lk_mpp ) THEN
635	zbathy(:,:) = FLOAT( mbathy(:,:) )
636	CALL lbc_lnk( zbathy, 'T', 1. )
637	mbathy(:,:) = INT( zbathy(:,:) )
638	ENDIF
639
640	! ! East-west cyclic boundary conditions
641	IF( nperio == 0 ) THEN
642	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
643	IF( lk_mpp ) THEN
644	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
645	IF( jperio /= 1 ) mbathy(1,:) = 0
646	ENDIF
647	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
648	IF( jperio /= 1 ) mbathy(nlci,:) = 0
649	ENDIF
650	ELSE
651	IF( ln_zco .OR. ln_zps ) THEN
652	mbathy( 1 ,:) = 0
653	mbathy(jpi,:) = 0
654	ELSE
655	mbathy( 1 ,:) = jpkm1
656	mbathy(jpi,:) = jpkm1
657	ENDIF
658	ENDIF
659	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
660	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
661	mbathy( 1 ,:) = mbathy(jpim1,:)
662	mbathy(jpi,:) = mbathy( 2 ,:)
663	ELSEIF( nperio == 2 ) THEN
664	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
665	ELSE
666	IF(lwp) WRITE(numout,*) ' e r r o r'
667	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
668	! STOP 'dom_mba'
669	ENDIF
670
671	! Boundary condition on mbathy
672	IF( .NOT.lk_mpp ) THEN
673	!!gm !!bug ??? think about it !
674	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
675	zbathy(:,:) = FLOAT( mbathy(:,:) )
676	CALL lbc_lnk( zbathy, 'T', 1. )
677	mbathy(:,:) = INT( zbathy(:,:) )
678	ENDIF
679
680	! Number of ocean level inferior or equal to jpkm1
681	ikmax = 0
682	DO jj = 1, jpj
683	DO ji = 1, jpi
684	ikmax = MAX( ikmax, mbathy(ji,jj) )
685	END DO
686	END DO
687	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
688	IF( ikmax > jpkm1 ) THEN
689	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
690	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
691	ELSE IF( ikmax < jpkm1 ) THEN
692	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
693	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
694	ENDIF
695
696	IF( lwp .AND. nprint == 1 ) THEN ! control print
697	WRITE(numout,*)
698	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
699	WRITE(numout,*) ' ------------------'
700	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
701	WRITE(numout,*)
702	ENDIF
703	!
704	END SUBROUTINE zgr_bat_ctl
705
706
707	SUBROUTINE zgr_zco
708	!!----------------------------------------------------------------------
709	!! * ROUTINE zgr_zco *
710	!!
711	!! ** Purpose : define the z-coordinate system
712	!!
713	!! ** Method : set 3D coord. arrays to reference 1D array
714	!!----------------------------------------------------------------------
715	INTEGER :: jk
716	!!----------------------------------------------------------------------
717	!
718	DO jk = 1, jpk
719	fsdept(:,:,jk) = gdept_0(jk)
720	fsdepw(:,:,jk) = gdepw_0(jk)
721	fsde3w(:,:,jk) = gdepw_0(jk)
722	fse3t (:,:,jk) = e3t_0(jk)
723	fse3u (:,:,jk) = e3t_0(jk)
724	fse3v (:,:,jk) = e3t_0(jk)
725	fse3f (:,:,jk) = e3t_0(jk)
726	fse3w (:,:,jk) = e3w_0(jk)
727	fse3uw(:,:,jk) = e3w_0(jk)
728	fse3vw(:,:,jk) = e3w_0(jk)
729	END DO
730	!
731	END SUBROUTINE zgr_zco
732
733	!!----------------------------------------------------------------------
734	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
735	!!----------------------------------------------------------------------
736
737	SUBROUTINE zgr_zps
738	!!----------------------------------------------------------------------
739	!! * ROUTINE zgr_zps *
740	!!
741	!! ** Purpose : the depth and vertical scale factor in partial step
742	!! z-coordinate case
743	!!
744	!! ** Method : Partial steps : computes the 3D vertical scale factors
745	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
746	!! a partial step representation of bottom topography.
747	!!
748	!! The reference depth of model levels is defined from an analytical
749	!! function the derivative of which gives the reference vertical
750	!! scale factors.
751	!! From depth and scale factors reference, we compute there new value
752	!! with partial steps on 3d arrays ( i, j, k ).
753	!!
754	!! w-level: gdepw(i,j,k) = fsdep(k)
755	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
756	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
757	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
758	!!
759	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
760	!! we find the mbathy index of the depth at each grid point.
761	!! This leads us to three cases:
762	!!
763	!! - bathy = 0 => mbathy = 0
764	!! - 1 < mbathy < jpkm1
765	!! - bathy > gdepw(jpk) => mbathy = jpkm1
766	!!
767	!! Then, for each case, we find the new depth at t- and w- levels
768	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
769	!! and f-points.
770	!!
771	!! This routine is given as an example, it must be modified
772	!! following the user s desiderata. nevertheless, the output as
773	!! well as the way to compute the model levels and scale factors
774	!! must be respected in order to insure second order accuracy
775	!! schemes.
776	!!
777	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
778	!! - - - - - - - gdept, gdepw and e3. are positives
779	!!
780	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
781	!!----------------------------------------------------------------------
782	INTEGER :: ji, jj, jk ! dummy loop indices
783	INTEGER :: ik, it ! temporary integers
784	INTEGER :: nlbelow ! temporary integer
785	LOGICAL :: ll_print ! Allow control print for debugging
786	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
787	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
788	REAL(wp) :: zmax, zmin ! Maximum and minimum depth
789	REAL(wp) :: zdiff ! temporary scalar
790	REAL(wp) :: zrefdep ! temporary scalar
791	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace
792	!!---------------------------------------------------------------------
793
794	IF(lwp) WRITE(numout,*)
795	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
796	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
797	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
798
799	ll_print=.FALSE. ! Local variable for debugging
800	!! ll_print=.TRUE.
801
802	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
803	WRITE(numout,*)
804	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
805	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
806	ENDIF
807
808
809	! bathymetry in level (from bathy_meter)
810	! ===================
811	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
812
813	zrefdep = 25.
814	nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below zrefdep
815	IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
816	zmin = gdepw_0(nlbelow) ! minimum depth = nlbelow-1 levels
817
818	mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available
819
820	! ! storage of land and island's number (zera and negative values) in mbathy
821	WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )
822
823	! bounded value of bathy
824	!!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0
825	DO jj = 1, jpj
826	DO ji= 1, jpi
827	IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0
828	ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )
829	ENDIF
830	END DO
831	END DO
832
833	! Compute mbathy for ocean points (i.e. the number of ocean levels)
834	! find the number of ocean levels such that the last level thickness
835	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
836	! e3t_0 is the reference level thickness
837	DO jk = jpkm1, 1, -1
838	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
839	DO jj = 1, jpj
840	DO ji = 1, jpi
841	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
842	END DO
843	END DO
844	END DO
845
846	! Scale factors and depth at T- and W-points
847	DO jk = 1, jpk ! intitialization to the reference z-coordinate
848	gdept(:,:,jk) = gdept_0(jk)
849	gdepw(:,:,jk) = gdepw_0(jk)
850	e3t (:,:,jk) = e3t_0 (jk)
851	e3w (:,:,jk) = e3w_0 (jk)
852	END DO
853	hdept(:,:) = gdept(:,:,2 )
854	hdepw(:,:) = gdepw(:,:,3 )
855	!
856	DO jj = 1, jpj
857	DO ji = 1, jpi
858	ik = mbathy(ji,jj)
859	IF( ik > 0 ) THEN ! ocean point only
860	! max ocean level case
861	IF( ik == jpkm1 ) THEN
862	zdepwp = bathy(ji,jj)
863	ze3tp = bathy(ji,jj) - gdepw_0(ik)
864	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
865	e3t(ji,jj,ik ) = ze3tp
866	e3t(ji,jj,ik+1) = ze3tp
867	e3w(ji,jj,ik ) = ze3wp
868	e3w(ji,jj,ik+1) = ze3tp
869	gdepw(ji,jj,ik+1) = zdepwp
870	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
871	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
872	!
873	ELSE ! standard case
874	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
875	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
876	ENDIF
877	!gm Bug? check the gdepw_0
878	! ... on ik
879	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
880	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
881	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
882	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
883	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
884	e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
885	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
886	! ... on ik+1
887	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
888	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
889	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
890	ENDIF
891	ENDIF
892	END DO
893	END DO
894	!
895	it = 0
896	DO jj = 1, jpj
897	DO ji = 1, jpi
898	ik = mbathy(ji,jj)
899	IF( ik > 0 ) THEN ! ocean point only
900	hdept(ji,jj) = gdept(ji,jj,ik )
901	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
902	e3tp (ji,jj) = e3t(ji,jj,ik )
903	e3wp (ji,jj) = e3w(ji,jj,ik )
904	! test
905	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
906	IF( zdiff <= 0. .AND. lwp ) THEN
907	it = it + 1
908	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
909	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
910	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
911	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
912	ENDIF
913	ENDIF
914	END DO
915	END DO
916
917	! Scale factors and depth at U-, V-, UW and VW-points
918	DO jk = 1, jpk ! initialisation to z-scale factors
919	e3u (:,:,jk) = e3t_0(jk)
920	e3v (:,:,jk) = e3t_0(jk)
921	e3uw(:,:,jk) = e3w_0(jk)
922	e3vw(:,:,jk) = e3w_0(jk)
923	END DO
924	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
925	DO jj = 1, jpjm1
926	DO ji = 1, fs_jpim1 ! vector opt.
927	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
928	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
929	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
930	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
931	END DO
932	END DO
933	END DO
934	! ! Boundary conditions
935	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
936	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
937	!
938	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
939	WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk)
940	WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk)
941	WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk)
942	WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk)
943	END DO
944
945	! Scale factor at F-point
946	DO jk = 1, jpk ! initialisation to z-scale factors
947	e3f (:,:,jk) = e3t_0(jk)
948	END DO
949	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
950	DO jj = 1, jpjm1
951	DO ji = 1, fs_jpim1 ! vector opt.
952	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
953	END DO
954	END DO
955	END DO
956	CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions
957	!
958	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
959	WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk)
960	END DO
961	!!gm bug ? : must be a do loop with mj0,mj1
962	!
963	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
964	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
965	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
966	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
967	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
968
969	! Control of the sign
970	IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
971	IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
972	IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
973	IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
974
975	! Compute gdep3w (vertical sum of e3w)
976	gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
977	DO jk = 2, jpk
978	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
979	END DO
980
981	! ! ================= !
982	IF(lwp .AND. ll_print) THEN ! Control print !
983	! ! ================= !
984	DO jj = 1,jpj
985	DO ji = 1, jpi
986	ik = MAX( mbathy(ji,jj), 1 )
987	zprt(ji,jj,1) = e3t (ji,jj,ik)
988	zprt(ji,jj,2) = e3w (ji,jj,ik)
989	zprt(ji,jj,3) = e3u (ji,jj,ik)
990	zprt(ji,jj,4) = e3v (ji,jj,ik)
991	zprt(ji,jj,5) = e3f (ji,jj,ik)
992	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
993	END DO
994	END DO
995	WRITE(numout,*)
996	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
997	WRITE(numout,*)
998	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
999	WRITE(numout,*)
1000	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1001	WRITE(numout,*)
1002	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1003	WRITE(numout,*)
1004	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1005	WRITE(numout,*)
1006	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1007	ENDIF
1008
1009	! ! =============== !
1010	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
1011	! ! =============== !
1012
1013	#if ! defined key_c1d
1014	! ! =================== !
1015	CALL zgr_bat_ctl ! Bathymetry check !
1016	! ! =================== !
1017	#endif
1018	END SUBROUTINE zgr_zps
1019
1020
1021	FUNCTION fssig( pk ) RESULT( pf )
1022	!!----------------------------------------------------------------------
1023	!! * ROUTINE eos_init *
1024	!!
1025	!! ** Purpose : provide the analytical function in s-coordinate
1026	!!
1027	!! ** Method : the function provide the non-dimensional position of
1028	!! T and W (i.e. between 0 and 1)
1029	!! T-points at integer values (between 1 and jpk)
1030	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1031	!!
1032	!! Reference : ???
1033	!!----------------------------------------------------------------------
1034	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
1035	REAL(wp) :: pf ! sigma value
1036	!!----------------------------------------------------------------------
1037	!
1038	pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) &
1039	& - TANH( rn_thetb * rn_theta ) ) &
1040	& * ( COSH( rn_theta ) &
1041	& + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) ) ) &
1042	& / ( 2.e0 * SINH( rn_theta ) )
1043	!
1044	END FUNCTION fssig
1045
1046
1047	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1048	!!----------------------------------------------------------------------
1049	!! * ROUTINE eos_init *
1050	!!
1051	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
1052	!!
1053	!! ** Method : the function provides the non-dimensional position of
1054	!! T and W (i.e. between 0 and 1)
1055	!! T-points at integer values (between 1 and jpk)
1056	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1057	!!
1058	!! Reference : ???
1059	!!----------------------------------------------------------------------
1060	REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate
1061	REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient
1062	REAL(wp) :: pf1 ! sigma value
1063	!!----------------------------------------------------------------------
1064	!
1065	IF ( rn_theta == 0 ) then ! uniform sigma
1066	pf1 = -(pk1-0.5) / REAL( jpkm1 )
1067	ELSE ! stretched sigma
1068	pf1 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + &
1069	& pbb * ( (tanh( rn_theta( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5rn_theta) ) / &
1070	& (2tanh(0.5rn_theta) ) )
1071	ENDIF
1072	!
1073	END FUNCTION fssig1
1074
1075
1076	SUBROUTINE zgr_sco
1077	!!----------------------------------------------------------------------
1078	!! * ROUTINE zgr_sco *
1079	!!
1080	!! ** Purpose : define the s-coordinate system
1081	!!
1082	!! ** Method : s-coordinate
1083	!! The depth of model levels is defined as the product of an
1084	!! analytical function by the local bathymetry, while the vertical
1085	!! scale factors are defined as the product of the first derivative
1086	!! of the analytical function by the bathymetry.
1087	!! (this solution save memory as depth and scale factors are not
1088	!! 3d fields)
1089	!! - Read bathymetry (in meters) at t-point and compute the
1090	!! bathymetry at u-, v-, and f-points.
1091	!! hbatu = mi( hbatt )
1092	!! hbatv = mj( hbatt )
1093	!! hbatf = mi( mj( hbatt ) )
1094	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
1095	!! function and its derivative given as function.
1096	!! gsigt(k) = fssig (k )
1097	!! gsigw(k) = fssig (k-0.5)
1098	!! esigt(k) = fsdsig(k )
1099	!! esigw(k) = fsdsig(k-0.5)
1100	!! This routine is given as an example, it must be modified
1101	!! following the user s desiderata. nevertheless, the output as
1102	!! well as the way to compute the model levels and scale factors
1103	!! must be respected in order to insure second order a!!uracy
1104	!! schemes.
1105	!!
1106	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1107	!!----------------------------------------------------------------------
1108	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1109	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1110	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars
1111	REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace
1112	REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - -
1113	!!
1114	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3
1115	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigt3
1116	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsi3w3
1117	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigt3
1118	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigw3
1119	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtu3
1120	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtv3
1121	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtf3
1122	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwu3
1123	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwv3
1124	!!
1125	NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
1126	!!----------------------------------------------------------------------
1127
1128	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
1129	READ ( numnam, namzgr_sco )
1130
1131	IF(lwp) THEN ! control print
1132	WRITE(numout,*)
1133	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1134	WRITE(numout,*) '~~~~~~~~~~~'
1135	WRITE(numout,*) ' Namelist namzgr_sco'
1136	WRITE(numout,*) ' sigma-stretching coeffs '
1137	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max
1138	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min
1139	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta
1140	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb
1141	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax
1142	WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma
1143	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb
1144	WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc
1145	ENDIF
1146
1147	gsigw3 = 0.0d0 ; gsigt3 = 0.0d0 ; gsi3w3 = 0.0d0
1148	esigt3 = 0.0d0 ; esigw3 = 0.0d0
1149	esigtu3 = 0.0d0 ; esigtv3 = 0.0d0 ; esigtf3 = 0.0d0
1150	esigwu3 = 0.0d0 ; esigwv3 = 0.0d0
1151
1152	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1153	hifu(:,:) = rn_sbot_min
1154	hifv(:,:) = rn_sbot_min
1155	hiff(:,:) = rn_sbot_min
1156
1157	! ! set maximum ocean depth
1158	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1159
1160	DO jj = 1, jpj
1161	DO ji = 1, jpi
1162	IF (bathy(ji,jj) .gt. 0.0) THEN
1163	bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1164	ENDIF
1165	END DO
1166	END DO
1167	! ! =============================
1168	! ! Define the envelop bathymetry (hbatt)
1169	! ! =============================
1170	! use r-value to create hybrid coordinates
1171	DO jj = 1, jpj
1172	DO ji = 1, jpi
1173	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
1174	END DO
1175	END DO
1176	!
1177	! Smooth the bathymetry (if required)
1178	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
1179	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1180	!
1181	jl = 0
1182	zrmax = 1.e0
1183	! ! ================ !
1184	DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
1185	! ! ================ !
1186	jl = jl + 1
1187	zrmax = 0.e0
1188	zmsk(:,:) = 0.e0
1189	DO jj = 1, nlcj
1190	DO ji = 1, nlci
1191	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1192	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1193	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1194	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1195	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1196	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
1197	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1.0
1198	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
1199	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1.0
1200	END DO
1201	END DO
1202	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1203	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1204	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
1205	DO jj = 1, nlcj
1206	DO ji = 1, nlci
1207	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1208	END DO
1209	END DO
1210	!
1211	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1212	!
1213	DO jj = 1, nlcj
1214	DO ji = 1, nlci
1215	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1216	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1217	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1218	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1219	IF( zmsk(ji,jj) == 1.0 ) THEN
1220	ztmp(ji,jj) = ( &
1221	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1222	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
1223	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1224	& ) / ( &
1225	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1226	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
1227	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1228	& )
1229	ENDIF
1230	END DO
1231	END DO
1232	!
1233	DO jj = 1, nlcj
1234	DO ji = 1, nlci
1235	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1236	END DO
1237	END DO
1238	!
1239	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
1240	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
1241	DO jj = 1, nlcj
1242	DO ji = 1, nlci
1243	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
1244	END DO
1245	END DO
1246	! ! ================ !
1247	END DO ! End loop !
1248	! ! ================ !
1249	!
1250	! ! envelop bathymetry saved in hbatt
1251	hbatt(:,:) = zenv(:,:)
1252	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
1253	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1254	DO jj = 1, jpj
1255	DO ji = 1, jpi
1256	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1257	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1258	END DO
1259	END DO
1260	ENDIF
1261	!
1262	IF(lwp) THEN ! Control print
1263	WRITE(numout,*)
1264	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1265	WRITE(numout,*)
1266	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1267	IF( nprint == 1 ) THEN
1268	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1269	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1270	ENDIF
1271	ENDIF
1272
1273	! ! ==============================
1274	! ! hbatu, hbatv, hbatf fields
1275	! ! ==============================
1276	IF(lwp) THEN
1277	WRITE(numout,*)
1278	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1279	ENDIF
1280	hbatu(:,:) = rn_sbot_min
1281	hbatv(:,:) = rn_sbot_min
1282	hbatf(:,:) = rn_sbot_min
1283	DO jj = 1, jpjm1
1284	DO ji = 1, jpim1 ! NO vector opt.
1285	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1286	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1287	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1288	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1289	END DO
1290	END DO
1291	!
1292	! Apply lateral boundary condition
1293	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1294	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
1295	DO jj = 1, jpj
1296	DO ji = 1, jpi
1297	IF( hbatu(ji,jj) == 0.e0 ) THEN
1298	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = rn_sbot_min
1299	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1300	ENDIF
1301	END DO
1302	END DO
1303	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
1304	DO jj = 1, jpj
1305	DO ji = 1, jpi
1306	IF( hbatv(ji,jj) == 0.e0 ) THEN
1307	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = rn_sbot_min
1308	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1309	ENDIF
1310	END DO
1311	END DO
1312	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
1313	DO jj = 1, jpj
1314	DO ji = 1, jpi
1315	IF( hbatf(ji,jj) == 0.e0 ) THEN
1316	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = rn_sbot_min
1317	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1318	ENDIF
1319	END DO
1320	END DO
1321
1322	!!bug: key_helsinki a verifer
1323	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1324	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1325	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1326	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1327
1328	IF( nprint == 1 .AND. lwp ) THEN
1329	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1330	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1331	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1332	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1333	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1334	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1335	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1336	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1337	ENDIF
1338	!! helsinki
1339
1340	! ! =======================
1341	! ! s-ccordinate fields (gdep., e3.)
1342	! ! =======================
1343	!
1344	! non-dimensional "sigma" for model level depth at w- and t-levels
1345
1346	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
1347	! ! below rn_hc, with uniform sigma in shallower waters
1348	DO ji = 1, jpi
1349	DO jj = 1, jpj
1350
1351	IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
1352	DO jk = 1, jpk
1353	gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1354	gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1355	END DO
1356	ELSE ! shallow water, uniform sigma
1357	DO jk = 1, jpk
1358	gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1359	gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1360	END DO
1361	ENDIF
1362	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
1363
1364
1365	DO jk = 1, jpkm1
1366	esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
1367	esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
1368	END DO
1369	esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ))
1370	esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
1371
1372	! Coefficients for vertical depth as the sum of e3w scale factors
1373	gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
1374	DO jk = 2, jpk
1375	gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
1376	END DO
1377
1378	DO jk = 1, jpk
1379	zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
1380	zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
1381	gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigt3(ji,jj,jk)+rn_hczcoeft)
1382	gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigw3(ji,jj,jk)+rn_hczcoefw)
1383	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsi3w3(ji,jj,jk)+rn_hczcoeft)
1384	END DO
1385
1386	ENDDO ! for all jj's
1387	ENDDO ! for all ji's
1388
1389
1390	DO ji=1,jpi
1391	DO jj=1,jpj
1392
1393	DO jk = 1, jpk
1394	esigtu3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) ) / &
1395	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1396	esigtv3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji,jj+1)esigt3(ji,jj+1,jk) ) / &
1397	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1398	esigtf3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) + &
1399	hbatt(ji,jj+1)esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)esigt3(ji+1,jj+1,jk) ) / &
1400	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1401	esigwu3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji+1,jj)esigw3(ji+1,jj,jk) ) / &
1402	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1403	esigwv3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji,jj+1)esigw3(ji,jj+1,jk) ) / &
1404	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1405
1406	e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1407	e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1408	e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1409	e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1410	!
1411	e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1412	e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1413	e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
1414	END DO
1415
1416	ENDDO
1417	ENDDO
1418
1419	ELSE ! not ln_s_sigma
1420
1421	DO jk = 1, jpk
1422	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1423	gsigt(jk) = -fssig( REAL(jk,wp) )
1424	END DO
1425	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1426	!
1427	! Coefficients for vertical scale factors at w-, t- levels
1428	!!gm bug : define it from analytical function, not like juste bellow....
1429	!!gm or betteroffer the 2 possibilities....
1430	DO jk = 1, jpkm1
1431	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
1432	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
1433	END DO
1434	esigw( 1 ) = 2.0_wp * (gsigt(1) - gsigw(1))
1435	esigt(jpk) = 2.0_wp * (gsigt(jpk) - gsigw(jpk))
1436
1437	!!gm original form
1438	!!org DO jk = 1, jpk
1439	!!org esigt(jk)=fsdsig( FLOAT(jk) )
1440	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
1441	!!org END DO
1442	!!gm
1443	!
1444	! Coefficients for vertical depth as the sum of e3w scale factors
1445	gsi3w(1) = 0.5 * esigw(1)
1446	DO jk = 2, jpk
1447	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
1448	END DO
1449	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1450	DO jk = 1, jpk
1451	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1452	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1453	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1454	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1455	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoeft)
1456	END DO
1457	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1458	DO jj = 1, jpj
1459	DO ji = 1, jpi
1460	DO jk = 1, jpk
1461	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1462	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1463	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1464	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1465	!
1466	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1467	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1468	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1469	END DO
1470	END DO
1471	END DO
1472
1473	ENDIF ! ln_s_sigma
1474
1475
1476	!
1477	!! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
1478
1479	fsdept(:,:,:) = gdept (:,:,:)
1480	fsdepw(:,:,:) = gdepw (:,:,:)
1481	fsde3w(:,:,:) = gdep3w(:,:,:)
1482	fse3t (:,:,:) = e3t (:,:,:)
1483	fse3u (:,:,:) = e3u (:,:,:)
1484	fse3v (:,:,:) = e3v (:,:,:)
1485	fse3f (:,:,:) = e3f (:,:,:)
1486	fse3w (:,:,:) = e3w (:,:,:)
1487	fse3uw(:,:,:) = e3uw (:,:,:)
1488	fse3vw(:,:,:) = e3vw (:,:,:)
1489	!!
1490	! HYBRID :
1491	DO jj = 1, jpj
1492	DO ji = 1, jpi
1493	DO jk = 1, jpkm1
1494	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1495	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1496	END DO
1497	END DO
1498	END DO
1499	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1500	& ' MAX ', MAXVAL( mbathy(:,:) )
1501
1502
1503	! ! ===========
1504	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
1505	! ! ===========
1506
1507	#if ! defined key_c1d
1508	! ! =================== !
1509	CALL zgr_bat_ctl ! Bathymetry check !
1510	! ! =================== !
1511	#endif
1512
1513	! ! =============
1514	IF(lwp) THEN ! Control print
1515	! ! =============
1516	WRITE(numout,*)
1517	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1518	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
1519	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
1520	ENDIF
1521	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1522	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1523	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
1524	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
1525	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
1526	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
1527	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
1528	& ' w ', MINVAL( fse3w (:,:,:) )
1529
1530	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
1531	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
1532	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
1533	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
1534	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
1535	& ' w ', MAXVAL( fse3w (:,:,:) )
1536	ENDIF
1537	!
1538	IF(lwp) THEN ! selected vertical profiles
1539	WRITE(numout,*)
1540	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1541	WRITE(numout,*) ' ~~~~~~ --------------------'
1542	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1543	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
1544	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
1545	DO jj = mj0(20), mj1(20)
1546	DO ji = mi0(20), mi1(20)
1547	WRITE(numout,*)
1548	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1549	WRITE(numout,*) ' ~~~~~~ --------------------'
1550	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1551	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1552	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1553	END DO
1554	END DO
1555	DO jj = mj0(74), mj1(74)
1556	DO ji = mi0(100), mi1(100)
1557	WRITE(numout,*)
1558	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1559	WRITE(numout,*) ' ~~~~~~ --------------------'
1560	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
1561	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
1562	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
1563	END DO
1564	END DO
1565	ENDIF
1566
1567	!!gm bug? no more necessary? if ! defined key_helsinki
1568	DO jk = 1, jpk
1569	DO jj = 1, jpj
1570	DO ji = 1, jpi
1571	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1572	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1573	CALL ctl_stop( ctmp1 )
1574	ENDIF
1575	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1576	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1577	CALL ctl_stop( ctmp1 )
1578	ENDIF
1579	END DO
1580	END DO
1581	END DO
1582	!!gm bug #endif
1583	!
1584	END SUBROUTINE zgr_sco
1585
1586
1587	!!======================================================================
1588	END MODULE domzgr

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

Download in other formats: