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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4147

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

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

Download in other formats: