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

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

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

Context Navigation

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

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

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

Download in other formats: