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

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

domzgr.F90 in branches/dev_004_VVL/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/dev_004_VVL/NEMO/OPA_SRC/DOM/domzgr.F90 @ 1385

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

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

Download in other formats: