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/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9746

Last change on this file since 9746 was 9746, checked in by mathiot, 7 years ago
isf geometry: add comments and minor cleaning
File size: 125.8 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	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
20	!!----------------------------------------------------------------------
21
22	!!----------------------------------------------------------------------
23	!! dom_zgr : defined the ocean vertical coordinate system
24	!! zgr_bat : bathymetry fields (levels and meters)
25	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
26	!! zgr_bat_ctl : check the bathymetry files
27	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
28	!! zgr_z : reference z-coordinate
29	!! zgr_zco : z-coordinate
30	!! zgr_zps : z-coordinate with partial steps
31	!! zgr_sco : s-coordinate
32	!! fssig : tanh stretch function
33	!! fssig1 : Song and Haidvogel 1994 stretch function
34	!! fgamma : Siddorn and Furner 2012 stretching function
35	!!---------------------------------------------------------------------
36	USE oce ! ocean variables
37	USE dom_oce ! ocean domain
38	USE closea ! closed seas
39	USE c1d ! 1D vertical configuration
40	USE in_out_manager ! I/O manager
41	USE iom ! I/O library
42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
43	USE lib_mpp ! distributed memory computing library
44	USE wrk_nemo ! Memory allocation
45	USE timing ! Timing
46
47	IMPLICIT NONE
48	PRIVATE
49
50	PUBLIC dom_zgr ! called by dom_init.F90
51
52	! remove subglacial lake under the ice sheet.
53	LOGICAL , PUBLIC :: ln_isfsubgl
54	REAL(wp), PUBLIC :: rn_isfsubgllon, rn_isfsubgllat
55	! !!* Namelist namzgr_sco *
56	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
57	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
58	!
59	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
60	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
61	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
62	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
63	! Song and Haidvogel 1994 stretching parameters
64	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
65	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
66	REAL(wp) :: rn_bb ! stretching parameter
67	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
68	! Siddorn and Furner stretching parameters
69	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
70	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
71	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
72	REAL(wp) :: rn_zs ! depth of surface grid box
73	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
74	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
75	REAL(wp) :: rn_zb_b ! offset for calculating Zb
76
77	!! * Substitutions
78	# include "domzgr_substitute.h90"
79	# include "vectopt_loop_substitute.h90"
80	!!----------------------------------------------------------------------
81	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
82	!! $Id$
83	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
84	!!----------------------------------------------------------------------
85	CONTAINS
86
87	SUBROUTINE dom_zgr
88	!!----------------------------------------------------------------------
89	!! * ROUTINE dom_zgr *
90	!!
91	!! ** Purpose : set the depth of model levels and the resulting
92	!! vertical scale factors.
93	!!
94	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
95	!! - read/set ocean depth and ocean levels (bathy, mbathy)
96	!! - vertical coordinate (gdep., e3.) depending on the
97	!! coordinate chosen :
98	!! ln_zco=T z-coordinate
99	!! ln_zps=T z-coordinate with partial steps
100	!! ln_zco=T s-coordinate
101	!!
102	!! ** Action : define gdep., e3., mbathy and bathy
103	!!----------------------------------------------------------------------
104	INTEGER :: ioptio, ibat ! local integer
105	INTEGER :: ios
106	!
107	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
108	!!----------------------------------------------------------------------
109	!
110	IF( nn_timing == 1 ) CALL timing_start('dom_zgr')
111	!
112	REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
113	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
114	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
115
116	REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
117	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
118	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
119	IF(lwm) WRITE ( numond, namzgr )
120
121	IF(lwp) THEN ! Control print
122	WRITE(numout,*)
123	WRITE(numout,*) 'dom_zgr : vertical coordinate'
124	WRITE(numout,*) '~~~~~~~'
125	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
126	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
127	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
128	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
129	WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
130	ENDIF
131
132	ioptio = 0 ! Check Vertical coordinate options
133	IF( ln_zco ) ioptio = ioptio + 1
134	IF( ln_zps ) ioptio = ioptio + 1
135	IF( ln_sco ) ioptio = ioptio + 1
136	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
137	!
138	! Build the vertical coordinate system
139	! ------------------------------------
140	CALL zgr_z ! Reference z-coordinate system (always called)
141	CALL zgr_bat ! Bathymetry fields (levels and meters)
142	IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain
143	IF( ln_zco ) CALL zgr_zco ! z-coordinate
144	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
145	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
146	!
147	! final adjustment of mbathy & check
148	! -----------------------------------
149	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
150	IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
151	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
152	CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
153	!
154	IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain
155	ibat = mbathy(2,2)
156	mbathy(:,:) = ibat
157	END IF
158	!
159	IF( nprint == 1 .AND. lwp ) THEN
160	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
161	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
162	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
163	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
164	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
165	& ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)), &
166	& ' w ', MINVAL( e3w_0(:,:,:) )
167
168	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
169	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
170	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
171	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
172	& ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)), &
173	& ' w ', MAXVAL( e3w_0(:,:,:) )
174	ENDIF
175	!
176	IF( nn_timing == 1 ) CALL timing_stop('dom_zgr')
177	!
178	END SUBROUTINE dom_zgr
179
180
181	SUBROUTINE zgr_z
182	!!----------------------------------------------------------------------
183	!! * ROUTINE zgr_z *
184	!!
185	!! ** Purpose : set the depth of model levels and the resulting
186	!! vertical scale factors.
187	!!
188	!! ** Method : z-coordinate system (use in all type of coordinate)
189	!! The depth of model levels is defined from an analytical
190	!! function the derivative of which gives the scale factors.
191	!! both depth and scale factors only depend on k (1d arrays).
192	!! w-level: gdepw_1d = gdep(k)
193	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
194	!! t-level: gdept_1d = gdep(k+0.5)
195	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
196	!!
197	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
198	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
199	!!
200	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
201	!!----------------------------------------------------------------------
202	INTEGER :: jk ! dummy loop indices
203	REAL(wp) :: zt, zw ! temporary scalars
204	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
205	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
206	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
207	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
208	!!----------------------------------------------------------------------
209	!
210	IF( nn_timing == 1 ) CALL timing_start('zgr_z')
211	!
212	! Set variables from parameters
213	! ------------------------------
214	zkth = ppkth ; zacr = ppacr
215	zdzmin = ppdzmin ; zhmax = pphmax
216	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
217
218	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
219	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
220	IF( ppa1 == pp_to_be_computed .AND. &
221	& ppa0 == pp_to_be_computed .AND. &
222	& ppsur == pp_to_be_computed ) THEN
223	!
224	#if defined key_agrif
225	za1 = ( ppdzmin - pphmax / FLOAT(jpkdta-1) ) &
226	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * ( LOG( COSH( (jpkdta - ppkth) / ppacr) )&
227	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
228	#else
229	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
230	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
231	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
232	#endif
233	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
234	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
235	ELSE
236	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
237	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
238	ENDIF
239
240	IF(lwp) THEN ! Parameter print
241	WRITE(numout,*)
242	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
243	WRITE(numout,*) ' ~~~~~~~'
244	IF( ppkth == 0._wp ) THEN
245	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
246	WRITE(numout,*) ' Total depth :', zhmax
247	#if defined key_agrif
248	WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1)
249	#else
250	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
251	#endif
252	ELSE
253	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
254	WRITE(numout,*) ' zsur, za0, za1 computed from '
255	WRITE(numout,*) ' zdzmin = ', zdzmin
256	WRITE(numout,*) ' zhmax = ', zhmax
257	ENDIF
258	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
259	WRITE(numout,*) ' zsur = ', zsur
260	WRITE(numout,*) ' za0 = ', za0
261	WRITE(numout,*) ' za1 = ', za1
262	WRITE(numout,*) ' zkth = ', zkth
263	WRITE(numout,*) ' zacr = ', zacr
264	IF( ldbletanh ) THEN
265	WRITE(numout,*) ' (Double tanh za2 = ', za2
266	WRITE(numout,*) ' parameters) zkth2= ', zkth2
267	WRITE(numout,*) ' zacr2= ', zacr2
268	ENDIF
269	ENDIF
270	ENDIF
271
272
273	! Reference z-coordinate (depth - scale factor at T- and W-points)
274	! ======================
275	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
276	#if defined key_agrif
277	za1 = zhmax / FLOAT(jpkdta-1)
278	#else
279	za1 = zhmax / FLOAT(jpk-1)
280	#endif
281	DO jk = 1, jpk
282	zw = FLOAT( jk )
283	zt = FLOAT( jk ) + 0.5_wp
284	gdepw_1d(jk) = ( zw - 1 ) * za1
285	gdept_1d(jk) = ( zt - 1 ) * za1
286	e3w_1d (jk) = za1
287	e3t_1d (jk) = za1
288	END DO
289	ELSE ! Madec & Imbard 1996 function
290	IF( .NOT. ldbletanh ) THEN
291	DO jk = 1, jpk
292	zw = REAL( jk , wp )
293	zt = REAL( jk , wp ) + 0.5_wp
294	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
295	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
296	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
297	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
298	END DO
299	ELSE
300	DO jk = 1, jpk
301	zw = FLOAT( jk )
302	zt = FLOAT( jk ) + 0.5_wp
303	! Double tanh function
304	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
305	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
306	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
307	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
308	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
309	& + za2 * TANH( (zw-zkth2) / zacr2 )
310	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
311	& + za2 * TANH( (zt-zkth2) / zacr2 )
312	END DO
313	ENDIF
314	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
315	ENDIF
316
317	IF ( ln_isfcav ) THEN
318	! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
319	! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
320	DO jk = 1, jpkm1
321	e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
322	END DO
323	e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
324
325	DO jk = 2, jpk
326	e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
327	END DO
328	e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
329	END IF
330
331	!!gm BUG in s-coordinate this does not work!
332	! deepest/shallowest W level Above/Below ~10m
333	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
334	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
335	nla10 = nlb10 - 1 ! deepest W level Above ~10m
336	!!gm end bug
337
338	IF(lwp) THEN ! control print
339	WRITE(numout,*)
340	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
341	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
342	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
343	ENDIF
344	DO jk = 1, jpk ! control positivity
345	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
346	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
347	END DO
348	!
349	IF( nn_timing == 1 ) CALL timing_stop('zgr_z')
350	!
351	END SUBROUTINE zgr_z
352
353
354	SUBROUTINE zgr_bat
355	!!----------------------------------------------------------------------
356	!! * ROUTINE zgr_bat *
357	!!
358	!! ** Purpose : set bathymetry both in levels and meters
359	!!
360	!! ** Method : read or define mbathy and bathy arrays
361	!! * level bathymetry:
362	!! The ocean basin geometry is given by a two-dimensional array,
363	!! mbathy, which is defined as follow :
364	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
365	!! at t-point (ji,jj).
366	!! = 0 over the continental t-point.
367	!! The array mbathy is checked to verified its consistency with
368	!! model option. in particular:
369	!! mbathy must have at least 1 land grid-points (mbathy<=0)
370	!! along closed boundary.
371	!! mbathy must be cyclic IF jperio=1.
372	!! mbathy must be lower or equal to jpk-1.
373	!! isolated ocean grid points are suppressed from mbathy
374	!! since they are only connected to remaining
375	!! ocean through vertical diffusion.
376	!! ntopo=-1 : rectangular channel or bassin with a bump
377	!! ntopo= 0 : flat rectangular channel or basin
378	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
379	!! bathy is read in 'bathy_meter.nc' NetCDF file
380	!!
381	!! ** Action : - mbathy: level bathymetry (in level index)
382	!! - bathy : meter bathymetry (in meters)
383	!!----------------------------------------------------------------------
384	INTEGER :: ji, jj, jl, jk ! dummy loop indices
385	INTEGER :: inum ! temporary logical unit
386	INTEGER :: ierror ! error flag
387	INTEGER :: ii_bump, ij_bump, ih ! bump center position
388	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
389	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
390	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
391	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data
392	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data
393	!!----------------------------------------------------------------------
394	!
395	IF( nn_timing == 1 ) CALL timing_start('zgr_bat')
396	!
397	IF(lwp) WRITE(numout,*)
398	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
399	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
400	!
401	! (ISF) initialisation ice shelf draft and top level
402	risfdep(:,:)=0._wp
403	misfdep(:,:)=1
404	! ! ================== !
405	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
406	! ! ================== !
407	! ! global domain level and meter bathymetry (idta,zdta)
408	!
409	ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
410	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
411	ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
412	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
413	!
414	IF( ntopo == 0 ) THEN ! flat basin
415	IF(lwp) WRITE(numout,*)
416	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
417	IF( rn_bathy > 0.01 ) THEN
418	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
419	zdta(:,:) = rn_bathy
420	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
421	idta(:,:) = jpkm1
422	ELSE ! z-coordinate (zco or zps): step-like topography
423	idta(:,:) = jpkm1
424	DO jk = 1, jpkm1
425	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
426	END DO
427	ENDIF
428	ELSE
429	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
430	idta(:,:) = jpkm1 ! before last level
431	zdta(:,:) = gdepw_1d(jpk) ! last w-point depth
432	h_oce = gdepw_1d(jpk)
433	ENDIF
434	ELSE ! bump centered in the basin
435	IF(lwp) WRITE(numout,*)
436	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
437	ii_bump = jpidta / 2 ! i-index of the bump center
438	ij_bump = jpjdta / 2 ! j-index of the bump center
439	r_bump = 50000._wp ! bump radius (meters)
440	h_bump = 2700._wp ! bump height (meters)
441	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
442	IF(lwp) WRITE(numout,*) ' bump characteristics: '
443	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
444	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
445	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
446	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
447	!
448	DO jj = 1, jpjdta ! zdta :
449	DO ji = 1, jpidta
450	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
451	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
452	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
453	END DO
454	END DO
455	! ! idta :
456	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
457	idta(:,:) = jpkm1
458	ELSE ! z-coordinate (zco or zps): step-like topography
459	idta(:,:) = jpkm1
460	DO jk = 1, jpkm1
461	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
462	END DO
463	ENDIF
464	ENDIF
465	! ! set GLOBAL boundary conditions
466	! ! Caution : idta on the global domain: use of jperio, not nperio
467	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
468	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
469	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
470	ELSEIF( jperio == 2 ) THEN
471	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
472	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
473	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
474	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp
475	ELSE
476	ih = 0 ; zh = 0._wp
477	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
478	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
479	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
480	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
481	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
482	ENDIF
483
484	! ! local domain level and meter bathymetries (mbathy,bathy)
485	mbathy(:,:) = 0 ! set to zero extra halo points
486	bathy (:,:) = 0._wp ! (require for mpp case)
487	DO jj = 1, nlcj ! interior values
488	DO ji = 1, nlci
489	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
490	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
491	END DO
492	END DO
493	!
494	DEALLOCATE( idta, zdta )
495	!
496	! ! ================ !
497	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
498	! ! ================ !
499	!
500	IF( ln_zco ) THEN ! zco : read level bathymetry
501	CALL iom_open ( 'bathy_level.nc', inum )
502	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
503	CALL iom_close( inum )
504	mbathy(:,:) = INT( bathy(:,:) )
505	! ! =====================
506	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
507	! ! =====================
508	IF( nn_cla == 0 ) THEN
509	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
510	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
511	DO ji = mi0(ii0), mi1(ii1)
512	DO jj = mj0(ij0), mj1(ij1)
513	mbathy(ji,jj) = 15
514	END DO
515	END DO
516	IF(lwp) WRITE(numout,*)
517	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
518	!
519	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
520	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
521	DO ji = mi0(ii0), mi1(ii1)
522	DO jj = mj0(ij0), mj1(ij1)
523	mbathy(ji,jj) = 12
524	END DO
525	END DO
526	IF(lwp) WRITE(numout,*)
527	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
528	ENDIF
529	!
530	ENDIF
531	!
532	ENDIF
533	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
534	CALL iom_open ( 'bathy_meter.nc', inum )
535	IF ( ln_isfcav ) THEN
536	CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
537	ELSE
538	CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr )
539	END IF
540	CALL iom_close( inum )
541	!
542	IF ( ln_isfcav ) THEN
543	CALL iom_open ( 'isf_draft_meter.nc', inum )
544	CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep )
545	CALL iom_close( inum )
546	WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
547	END IF
548	!
549	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
550	!
551	IF( nn_cla == 0 ) THEN
552	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
553	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
554	DO ji = mi0(ii0), mi1(ii1)
555	DO jj = mj0(ij0), mj1(ij1)
556	bathy(ji,jj) = 284._wp
557	END DO
558	END DO
559	IF(lwp) WRITE(numout,*)
560	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
561	!
562	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
563	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
564	DO ji = mi0(ii0), mi1(ii1)
565	DO jj = mj0(ij0), mj1(ij1)
566	bathy(ji,jj) = 137._wp
567	END DO
568	END DO
569	IF(lwp) WRITE(numout,*)
570	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
571	ENDIF
572	!
573	ENDIF
574	!
575	ENDIF
576	! ! =============== !
577	ELSE ! error !
578	! ! =============== !
579	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
580	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
581	ENDIF
582	!
583	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
584	!
585	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
586	! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
587	IF ( ln_isfcav ) THEN
588	WHERE (bathy == risfdep)
589	bathy = 0.0_wp ; risfdep = 0.0_wp
590	END WHERE
591	END IF
592	! end patch
593	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
594	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
595	ENDIF
596	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
597	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
598	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
599	END WHERE
600	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
601	ENDIF
602	!
603	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')
604	!
605	END SUBROUTINE zgr_bat
606
607
608	SUBROUTINE zgr_bat_zoom
609	!!----------------------------------------------------------------------
610	!! * ROUTINE zgr_bat_zoom *
611	!!
612	!! ** Purpose : - Close zoom domain boundary if necessary
613	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
614	!!
615	!! ** Method :
616	!!
617	!! ** Action : - update mbathy: level bathymetry (in level index)
618	!!----------------------------------------------------------------------
619	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
620	!!----------------------------------------------------------------------
621	!
622	IF(lwp) WRITE(numout,*)
623	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
624	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
625	!
626	! Zoom domain
627	! ===========
628	!
629	! Forced closed boundary if required
630	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
631	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
632	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
633	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
634	!
635	! Configuration specific domain modifications
636	! (here, ORCA arctic configuration: suppress Med Sea)
637	IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
638	SELECT CASE ( jp_cfg )
639	! ! =======================
640	CASE ( 2 ) ! ORCA_R2 configuration
641	! ! =======================
642	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
643	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
644	ij0 = 98 ; ij1 = 110
645	! ! =======================
646	CASE ( 05 ) ! ORCA_R05 configuration
647	! ! =======================
648	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
649	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
650	ij0 = 314 ; ij1 = 370
651	END SELECT
652	!
653	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
654	!
655	ENDIF
656	!
657	END SUBROUTINE zgr_bat_zoom
658
659
660	SUBROUTINE zgr_bat_ctl
661	!!----------------------------------------------------------------------
662	!! * ROUTINE zgr_bat_ctl *
663	!!
664	!! ** Purpose : check the bathymetry in levels
665	!!
666	!! ** Method : The array mbathy is checked to verified its consistency
667	!! with the model options. in particular:
668	!! mbathy must have at least 1 land grid-points (mbathy<=0)
669	!! along closed boundary.
670	!! mbathy must be cyclic IF jperio=1.
671	!! mbathy must be lower or equal to jpk-1.
672	!! isolated ocean grid points are suppressed from mbathy
673	!! since they are only connected to remaining
674	!! ocean through vertical diffusion.
675	!! C A U T I O N : mbathy will be modified during the initializa-
676	!! tion phase to become the number of non-zero w-levels of a water
677	!! column, with a minimum value of 1.
678	!!
679	!! ** Action : - update mbathy: level bathymetry (in level index)
680	!! - update bathy : meter bathymetry (in meters)
681	!!----------------------------------------------------------------------
682	!!
683	INTEGER :: ji, jj, jl ! dummy loop indices
684	INTEGER :: icompt, ibtest, ikmax ! temporary integers
685	REAL(wp), POINTER, DIMENSION(:,:) :: zbathy
686
687	!!----------------------------------------------------------------------
688	!
689	IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl')
690	!
691	CALL wrk_alloc( jpi, jpj, zbathy )
692	!
693	IF(lwp) WRITE(numout,*)
694	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
695	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
696	! ! Suppress isolated ocean grid points
697	IF(lwp) WRITE(numout,*)
698	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
699	IF(lwp) WRITE(numout,*)' -----------------------------------'
700	icompt = 0
701	DO jl = 1, 2
702	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
703	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
704	mbathy(jpi,:) = mbathy( 2 ,:)
705	ENDIF
706	DO jj = 2, jpjm1
707	DO ji = 2, jpim1
708	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
709	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
710	IF( ibtest < mbathy(ji,jj) ) THEN
711	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
712	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
713	mbathy(ji,jj) = ibtest
714	icompt = icompt + 1
715	ENDIF
716	END DO
717	END DO
718	END DO
719	IF( lk_mpp ) CALL mpp_sum( icompt )
720	IF( icompt == 0 ) THEN
721	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
722	ELSE
723	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
724	ENDIF
725	IF( lk_mpp ) THEN
726	zbathy(:,:) = FLOAT( mbathy(:,:) )
727	CALL lbc_lnk( zbathy, 'T', 1._wp )
728	mbathy(:,:) = INT( zbathy(:,:) )
729	ENDIF
730	! ! East-west cyclic boundary conditions
731	IF( nperio == 0 ) THEN
732	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
733	IF( lk_mpp ) THEN
734	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
735	IF( jperio /= 1 ) mbathy(1,:) = 0
736	ENDIF
737	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
738	IF( jperio /= 1 ) mbathy(nlci,:) = 0
739	ENDIF
740	ELSE
741	IF( ln_zco .OR. ln_zps ) THEN
742	mbathy( 1 ,:) = 0
743	mbathy(jpi,:) = 0
744	ELSE
745	mbathy( 1 ,:) = jpkm1
746	mbathy(jpi,:) = jpkm1
747	ENDIF
748	ENDIF
749	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
750	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
751	mbathy( 1 ,:) = mbathy(jpim1,:)
752	mbathy(jpi,:) = mbathy( 2 ,:)
753	ELSEIF( nperio == 2 ) THEN
754	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
755	ELSE
756	IF(lwp) WRITE(numout,*) ' e r r o r'
757	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
758	! STOP 'dom_mba'
759	ENDIF
760	! Boundary condition on mbathy
761	IF( .NOT.lk_mpp ) THEN
762	!!gm !!bug ??? think about it !
763	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
764	zbathy(:,:) = FLOAT( mbathy(:,:) )
765	CALL lbc_lnk( zbathy, 'T', 1._wp )
766	mbathy(:,:) = INT( zbathy(:,:) )
767	ENDIF
768	! Number of ocean level inferior or equal to jpkm1
769	ikmax = 0
770	DO jj = 1, jpj
771	DO ji = 1, jpi
772	ikmax = MAX( ikmax, mbathy(ji,jj) )
773	END DO
774	END DO
775	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
776	IF( ikmax > jpkm1 ) THEN
777	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
778	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
779	ELSE IF( ikmax < jpkm1 ) THEN
780	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
781	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
782	ENDIF
783
784	IF( lwp .AND. nprint == 1 ) THEN ! control print
785	WRITE(numout,*)
786	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
787	WRITE(numout,*) ' ------------------'
788	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
789	WRITE(numout,*)
790	ENDIF
791	!
792	CALL wrk_dealloc( jpi, jpj, zbathy )
793	!
794	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl')
795	!
796	END SUBROUTINE zgr_bat_ctl
797
798
799	SUBROUTINE zgr_bot_level
800	!!----------------------------------------------------------------------
801	!! * ROUTINE zgr_bot_level *
802	!!
803	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
804	!!
805	!! ** Method : computes from mbathy with a minimum value of 1 over land
806	!!
807	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
808	!! ocean level at t-, u- & v-points
809	!! (min value = 1 over land)
810	!!----------------------------------------------------------------------
811	!!
812	INTEGER :: ji, jj ! dummy loop indices
813	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
814	!!----------------------------------------------------------------------
815	!
816	IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level')
817	!
818	CALL wrk_alloc( jpi, jpj, zmbk )
819	!
820	IF(lwp) WRITE(numout,*)
821	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
822	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
823	!
824	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
825
826	! ! bottom k-index of W-level = mbkt+1
827	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
828	DO ji = 1, jpim1
829	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
830	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
831	END DO
832	END DO
833	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
834	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
835	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
836	!
837	CALL wrk_dealloc( jpi, jpj, zmbk )
838	!
839	IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level')
840	!
841	END SUBROUTINE zgr_bot_level
842
843	SUBROUTINE zgr_top_level
844	!!----------------------------------------------------------------------
845	!! * ROUTINE zgr_bot_level *
846	!!
847	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
848	!!
849	!! ** Method : computes from misfdep with a minimum value of 1
850	!!
851	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
852	!! ocean level at t-, u- & v-points
853	!! (min value = 1)
854	!!----------------------------------------------------------------------
855	!!
856	INTEGER :: ji, jj ! dummy loop indices
857	REAL(wp), POINTER, DIMENSION(:,:) :: zmik
858	!!----------------------------------------------------------------------
859	!
860	IF( nn_timing == 1 ) CALL timing_start('zgr_top_level')
861	!
862	CALL wrk_alloc( jpi, jpj, zmik )
863	!
864	IF(lwp) WRITE(numout,*)
865	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
866	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
867	!
868	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
869	! ! top k-index of W-level (=mikt)
870	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
871	DO ji = 1, jpim1
872	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
873	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
874	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
875	END DO
876	END DO
877
878	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
879	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
880	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
881	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
882	!
883	CALL wrk_dealloc( jpi, jpj, zmik )
884	!
885	IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level')
886	!
887	END SUBROUTINE zgr_top_level
888
889	SUBROUTINE zgr_zco
890	!!----------------------------------------------------------------------
891	!! * ROUTINE zgr_zco *
892	!!
893	!! ** Purpose : define the z-coordinate system
894	!!
895	!! ** Method : set 3D coord. arrays to reference 1D array
896	!!----------------------------------------------------------------------
897	INTEGER :: jk
898	!!----------------------------------------------------------------------
899	!
900	IF( nn_timing == 1 ) CALL timing_start('zgr_zco')
901	!
902	DO jk = 1, jpk
903	gdept_0 (:,:,jk) = gdept_1d(jk)
904	gdepw_0 (:,:,jk) = gdepw_1d(jk)
905	gdep3w_0(:,:,jk) = gdepw_1d(jk)
906	e3t_0 (:,:,jk) = e3t_1d (jk)
907	e3u_0 (:,:,jk) = e3t_1d (jk)
908	e3v_0 (:,:,jk) = e3t_1d (jk)
909	e3f_0 (:,:,jk) = e3t_1d (jk)
910	e3w_0 (:,:,jk) = e3w_1d (jk)
911	e3uw_0 (:,:,jk) = e3w_1d (jk)
912	e3vw_0 (:,:,jk) = e3w_1d (jk)
913	END DO
914	!
915	IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')
916	!
917	END SUBROUTINE zgr_zco
918
919
920	SUBROUTINE zgr_zps
921	!!----------------------------------------------------------------------
922	!! * ROUTINE zgr_zps *
923	!!
924	!! ** Purpose : the depth and vertical scale factor in partial step
925	!! z-coordinate case
926	!!
927	!! ** Method : Partial steps : computes the 3D vertical scale factors
928	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
929	!! a partial step representation of bottom topography.
930	!!
931	!! The reference depth of model levels is defined from an analytical
932	!! function the derivative of which gives the reference vertical
933	!! scale factors.
934	!! From depth and scale factors reference, we compute there new value
935	!! with partial steps on 3d arrays ( i, j, k ).
936	!!
937	!! w-level: gdepw_0(i,j,k) = gdep(k)
938	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
939	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
940	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
941	!!
942	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
943	!! we find the mbathy index of the depth at each grid point.
944	!! This leads us to three cases:
945	!!
946	!! - bathy = 0 => mbathy = 0
947	!! - 1 < mbathy < jpkm1
948	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
949	!!
950	!! Then, for each case, we find the new depth at t- and w- levels
951	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
952	!! and f-points.
953	!!
954	!! This routine is given as an example, it must be modified
955	!! following the user s desiderata. nevertheless, the output as
956	!! well as the way to compute the model levels and scale factors
957	!! must be respected in order to insure second order accuracy
958	!! schemes.
959	!!
960	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
961	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
962	!!
963	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
964	!!----------------------------------------------------------------------
965	!!
966	INTEGER :: ji, jj, jk ! dummy loop indices
967	INTEGER :: ik, it, ikb, ikt ! temporary integers
968	LOGICAL :: ll_print ! Allow control print for debugging
969	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
970	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
971	REAL(wp) :: zmax ! Maximum depth
972	REAL(wp) :: zdiff ! temporary scalar
973	REAL(wp) :: zrefdep ! temporary scalar
974	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
975	!!---------------------------------------------------------------------
976	!
977	IF( nn_timing == 1 ) CALL timing_start('zgr_zps')
978	!
979	CALL wrk_alloc( jpi, jpj, jpk, zprt )
980	!
981	IF(lwp) WRITE(numout,*)
982	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
983	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
984	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
985
986	ll_print = .FALSE. ! Local variable for debugging
987
988	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
989	WRITE(numout,*)
990	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
991	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
992	ENDIF
993
994
995	! bathymetry in level (from bathy_meter)
996	! ===================
997	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
998	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
999	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
1000	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
1001	END WHERE
1002
1003	! Compute mbathy for ocean points (i.e. the number of ocean levels)
1004	! find the number of ocean levels such that the last level thickness
1005	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1006	! e3t_1d is the reference level thickness
1007	DO jk = jpkm1, 1, -1
1008	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1009	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
1010	END DO
1011
1012	IF ( ln_isfcav ) CALL zgr_isf
1013
1014	! Scale factors and depth at T- and W-points
1015	DO jk = 1, jpk ! intitialization to the reference z-coordinate
1016	gdept_0(:,:,jk) = gdept_1d(jk)
1017	gdepw_0(:,:,jk) = gdepw_1d(jk)
1018	e3t_0 (:,:,jk) = e3t_1d (jk)
1019	e3w_0 (:,:,jk) = e3w_1d (jk)
1020	END DO
1021	!
1022	DO jj = 1, jpj
1023	DO ji = 1, jpi
1024	ik = mbathy(ji,jj)
1025	IF( ik > 0 ) THEN ! ocean point only
1026	! max ocean level case
1027	IF( ik == jpkm1 ) THEN
1028	zdepwp = bathy(ji,jj)
1029	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1030	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1031	e3t_0(ji,jj,ik ) = ze3tp
1032	e3t_0(ji,jj,ik+1) = ze3tp
1033	e3w_0(ji,jj,ik ) = ze3wp
1034	e3w_0(ji,jj,ik+1) = ze3tp
1035	gdepw_0(ji,jj,ik+1) = zdepwp
1036	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1037	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1038	!
1039	ELSE ! standard case
1040	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1041	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1042	ENDIF
1043	!gm Bug? check the gdepw_1d
1044	! ... on ik
1045	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1046	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1047	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1048	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1049	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1050	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1051	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1052	! ... on ik+1
1053	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1054	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1055	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1056	ENDIF
1057	ENDIF
1058	END DO
1059	END DO
1060	!
1061	it = 0
1062	DO jj = 1, jpj
1063	DO ji = 1, jpi
1064	ik = mbathy(ji,jj)
1065	IF( ik > 0 ) THEN ! ocean point only
1066	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1067	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1068	! test
1069	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1070	IF( zdiff <= 0._wp .AND. lwp ) THEN
1071	it = it + 1
1072	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1073	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1074	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1075	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1076	ENDIF
1077	ENDIF
1078	END DO
1079	END DO
1080	!
1081	IF ( ln_isfcav ) THEN
1082	! (ISF) Definition of e3t, u, v, w for ISF case
1083	DO jj = 1, jpj
1084	DO ji = 1, jpi
1085	ik = misfdep(ji,jj)
1086	IF( ik > 1 ) THEN ! ice shelf point only
1087	IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
1088	gdepw_0(ji,jj,ik) = risfdep(ji,jj)
1089	!gm Bug? check the gdepw_0
1090	! ... on ik
1091	gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
1092	& * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
1093	& / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
1094	e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
1095	e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
1096
1097	IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
1098	e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
1099	ENDIF
1100	! ... on ik / ik-1
1101	e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))
1102	e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1103	! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
1104	gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
1105	ENDIF
1106	END DO
1107	END DO
1108	!
1109	it = 0
1110	DO jj = 1, jpj
1111	DO ji = 1, jpi
1112	ik = misfdep(ji,jj)
1113	IF( ik > 1 ) THEN ! ice shelf point only
1114	e3tp (ji,jj) = e3t_0(ji,jj,ik )
1115	e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
1116	! test
1117	zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
1118	IF( zdiff <= 0. .AND. lwp ) THEN
1119	it = it + 1
1120	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1121	WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
1122	WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1123	WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
1124	ENDIF
1125	ENDIF
1126	END DO
1127	END DO
1128	END IF
1129	! END (ISF)
1130
1131	! Scale factors and depth at U-, V-, UW and VW-points
1132	DO jk = 1, jpk ! initialisation to z-scale factors
1133	e3u_0 (:,:,jk) = e3t_1d(jk)
1134	e3v_0 (:,:,jk) = e3t_1d(jk)
1135	e3uw_0(:,:,jk) = e3w_1d(jk)
1136	e3vw_0(:,:,jk) = e3w_1d(jk)
1137	END DO
1138	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1139	DO jj = 1, jpjm1
1140	DO ji = 1, fs_jpim1 ! vector opt.
1141	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1142	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1143	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1144	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1145	END DO
1146	END DO
1147	END DO
1148	IF ( ln_isfcav ) THEN
1149	! (ISF) define e3uw (adapted for 2 cells in the water column)
1150	DO jj = 2, jpjm1
1151	DO ji = 2, fs_jpim1 ! vector opt.
1152	ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
1153	ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
1154	IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &
1155	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )
1156	ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
1157	ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
1158	IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &
1159	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )
1160	END DO
1161	END DO
1162	END IF
1163
1164	CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1165	CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1166	!
1167	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1168	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1169	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1170	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1171	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1172	END DO
1173
1174	! Scale factor at F-point
1175	DO jk = 1, jpk ! initialisation to z-scale factors
1176	e3f_0(:,:,jk) = e3t_1d(jk)
1177	END DO
1178	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1179	DO jj = 1, jpjm1
1180	DO ji = 1, fs_jpim1 ! vector opt.
1181	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1182	END DO
1183	END DO
1184	END DO
1185	CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1186	!
1187	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1188	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1189	END DO
1190	!!gm bug ? : must be a do loop with mj0,mj1
1191	!
1192	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1193	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1194	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1195	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1196	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1197
1198	! Control of the sign
1199	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1200	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1201	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1202	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1203
1204	! Compute gdep3w_0 (vertical sum of e3w)
1205	IF ( ln_isfcav ) THEN ! if cavity
1206	WHERE (misfdep == 0) misfdep = 1
1207	DO jj = 1,jpj
1208	DO ji = 1,jpi
1209	gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
1210	DO jk = 2, misfdep(ji,jj)
1211	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1212	END DO
1213	IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
1214	DO jk = misfdep(ji,jj) + 1, jpk
1215	gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1216	END DO
1217	END DO
1218	END DO
1219	ELSE ! no cavity
1220	gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1221	DO jk = 2, jpk
1222	gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1223	END DO
1224	END IF
1225	! ! ================= !
1226	IF(lwp .AND. ll_print) THEN ! Control print !
1227	! ! ================= !
1228	DO jj = 1,jpj
1229	DO ji = 1, jpi
1230	ik = MAX( mbathy(ji,jj), 1 )
1231	zprt(ji,jj,1) = e3t_0 (ji,jj,ik)
1232	zprt(ji,jj,2) = e3w_0 (ji,jj,ik)
1233	zprt(ji,jj,3) = e3u_0 (ji,jj,ik)
1234	zprt(ji,jj,4) = e3v_0 (ji,jj,ik)
1235	zprt(ji,jj,5) = e3f_0 (ji,jj,ik)
1236	zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1237	END DO
1238	END DO
1239	WRITE(numout,*)
1240	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1241	WRITE(numout,*)
1242	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1243	WRITE(numout,*)
1244	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1245	WRITE(numout,*)
1246	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1247	WRITE(numout,*)
1248	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1249	WRITE(numout,*)
1250	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1251	ENDIF
1252	!
1253	CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1254	!
1255	IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
1256	!
1257	END SUBROUTINE zgr_zps
1258
1259	SUBROUTINE zgr_isf
1260	!!----------------------------------------------------------------------
1261	!! * ROUTINE zgr_isf *
1262	!!
1263	!! ** Purpose : check the bathymetry in levels
1264	!!
1265	!! ** Method : The water column have to contained at least 2 cells
1266	!! Bathymetry and isfdraft are modified (dig/close) to respect
1267	!! this criterion.
1268	!! 0.0 = definition of minal ice shelf draft
1269	!! digging and filling base on depth criterion only
1270	!! 1.0 = set iceshelf to the minimum depth allowed
1271	!! 1.1 = ground ice shelf if water column less than X m
1272	!! 1.2 = ensure a minimum thickness for iceshelf cavity
1273	!! 1.3 = remove channels and single point 'bay'
1274	!! 1.4 = close single isolated point
1275	!! compute level
1276	!! 2.0 = compute level
1277	!! 2.1 = ensure misfdep is not below bathymetry after step 2.0
1278	!! digging and filling base on level criterion only
1279	!! 3.0 = dig to fit the 2 water level criterion (closed pool possible after this step)
1280	!! 3.1 = dig enough to ensure connectivity of all the cell beneath an ice shelf
1281	!! (most of the close pool should be remove after this step)
1282	!! 3.2 = fill chimney
1283	!!
1284	!! ** Action : - test compatibility between isfdraft and bathy
1285	!! - bathy and isfdraft are modified
1286	!!----------------------------------------------------------------------
1287	!!
1288	INTEGER :: ji, jj, jk, jn
1289	INTEGER :: ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
1290	INTEGER :: zmbathyij, zmbathyip1, zmbathyim1, zmbathyjp1, zmbathyjm1
1291	INTEGER , POINTER, DIMENSION(:,:) :: zmisfdep, zmbathy ! 2D workspace (ISH)
1292	!
1293	REAL(wp) :: zdepth, vmskjp1, vmskjm1, umskip1, umskim1, umskip1_r, umskim1_r, vmskjp1_r, vmskjm1_r
1294	REAL(wp) :: zisfdep_min
1295	REAL(wp), POINTER, DIMENSION(:,:) :: zdummy, zmask, zrisfdep ! 2D workspace (ISH)
1296	!
1297	!!-----------------
1298	! NAMELIST
1299	INTEGER :: ios
1300	INTEGER :: nn_kisfmax = 999. !
1301	REAL(wp) :: rn_isfdep_min = 10.0_wp ! ice shelf minimal thickness
1302	REAL(wp) :: rn_glhw_min = 1.0e-3 ! threshold on hw to define grounding line water
1303	REAL(wp) :: rn_isfhw_min = 1.0e-3 ! threshold on hw to define isf draft into the cavity
1304	REAL(wp) :: rn_isfshallow = 0.0_wp ! threshold to define shallow ice shelf cavity
1305	REAL(wp) :: rn_zisfmax = 6000.0_wp ! maximun meter of ice we are allowed to dig to assure connectivity
1306	LOGICAL :: ln_isfcheminey= .FALSE. , ln_isfconnect= .FALSE. , ln_isfchannel= .FALSE. ! remove cheminey, assure connectivity, remove channel in water column (on z space not level space)
1307	!!---------------------------------------------------------------------
1308	NAMELIST/namzgr_isf/nn_kisfmax, rn_zisfmax, &
1309	& rn_isfdep_min, rn_glhw_min, rn_isfhw_min, &
1310	& ln_isfcheminey, ln_isfconnect, ln_isfchannel, ln_isfsubgl, rn_isfsubgllon, rn_isfsubgllat
1311	!
1312	IF( nn_timing == 1 ) CALL timing_start('zgr_isf')
1313	!
1314	CALL wrk_alloc( jpi, jpj, zdummy, zmask, zrisfdep)
1315	CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
1316	!
1317	!
1318	REWIND( numnam_ref ) ! Namelist namzgr_isf in reference namelist : ice shelf geometry definition
1319	READ ( numnam_ref, namzgr_isf, IOSTAT = ios, ERR = 901)
1320	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in reference namelist', lwp )
1321
1322	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : ice shelf geometry definition
1323	READ ( numnam_cfg, namzgr_isf, IOSTAT = ios, ERR = 902 )
1324	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_isf in configuration namelist', lwp )
1325	IF(lwm) WRITE ( numond, namzgr_isf )
1326	!
1327	IF (lwp) THEN
1328	WRITE(numout,*) ' nn_kisfmax = ',nn_kisfmax !
1329	WRITE(numout,*) ' rn_zisfmax = ',rn_zisfmax !
1330	WRITE(numout,*) ' rn_isfdep_min = ',rn_isfdep_min !
1331	WRITE(numout,*) ' rn_isfhw_min = ',rn_isfhw_min !
1332	WRITE(numout,*) ' rn_glhw_min = ',rn_glhw_min !
1333	WRITE(numout,*) ' ln_isfcheminey = ',ln_isfcheminey !
1334	WRITE(numout,*) ' ln_isfconnect = ',ln_isfconnect !
1335	WRITE(numout,*) ' ln_isfchannel = ',ln_isfchannel !
1336	END IF
1337	!
1338	! 0.0 variable definition of minimal ice shelf draft allowed
1339	! an ice shelf draft have to be larger than the first level thickness
1340	zisfdep_min = MAX(rn_isfdep_min,e3t_1d(1))
1341
1342	! 1.0 set iceshelf to the minimum depth allowed
1343	WHERE(risfdep(:,:) > 0.0_wp .AND. risfdep(:,:) < zisfdep_min)
1344	risfdep(:,:)=zisfdep_min
1345	END WHERE
1346	!
1347	! 1.1 ground ice shelf if water column less than X m => set the grounding
1348	! line position
1349	WHERE( bathy(:,:) - risfdep(:,:) < rn_glhw_min )
1350	risfdep(:,:)=0.0 ; misfdep(:,:) = 1
1351	bathy (:,:)=0.0 ; mbathy (:,:) = 0
1352	END WHERE
1353	!
1354	! 1.2 ensure a minimum thickness for iceshelf cavity => avoid to negative
1355	! e3t if ssh + sum(e3t*tmask) < 0
1356	DO jj = 1, jpj
1357	DO ji = 1, jpi
1358	IF (bathy(ji,jj)-risfdep(ji,jj) < rn_isfhw_min) THEN
1359	risfdep(ji,jj) = bathy(ji,jj) - rn_isfhw_min
1360	! sanity check on risfdep (if < zisfdep_min) => we ground it
1361	IF (risfdep(ji,jj) < zisfdep_min) THEN
1362	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1363	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1364	END IF
1365	END IF
1366	END DO
1367	END DO
1368	!
1369	! 1.3 Remove channels and single point 'bay'.
1370	! channel could be created if connectivity is enforced later.
1371	IF (ln_isfchannel) THEN
1372	zmask(:,:)=1._wp
1373	WHERE (mbathy(:,:)==0.0)
1374	zmask(:,:)=0._wp
1375	END WHERE
1376	DO jj = 2, jpjm1
1377	DO ji = 2, jpim1
1378	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1379	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = connexion
1380	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = connexion
1381	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = connexion
1382	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = connexion
1383	! zonal channel and single bay
1384	IF ((umskip1+umskim1>=1) .AND. (vmskjp1+vmskjm1==0)) THEN
1385	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1386	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1387	END IF
1388	! meridionnal channel and single bay
1389	IF ((vmskjp1+vmskjm1>=1) .AND. (umskip1+umskim1==0)) THEN
1390	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1391	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1392	END IF
1393	END IF
1394	END DO
1395	END DO
1396	! ensure halo correct
1397	IF( lk_mpp ) THEN
1398	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1399	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1400	CALL lbc_lnk( risfdep, 'T', 1. )
1401	CALL lbc_lnk( bathy , 'T', 1. )
1402	ENDIF
1403	END IF
1404	!
1405	! 1.4 Closed single isolated point
1406	zmask(:,:)=1._wp
1407	WHERE (mbathy(:,:)==0.0)
1408	zmask(:,:)=0._wp
1409	END WHERE
1410	DO jj = 2, jpjm1
1411	DO ji = 2, jpim1
1412	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1413	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = connexion
1414	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = connexion
1415	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = connexion
1416	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = connexion
1417	! remove single point
1418	IF (umskip1+umskim1+vmskjp1+vmskjm1==0.0_wp) THEN
1419	risfdep(ji,jj)=0.0 ; misfdep(ji,jj) = 1
1420	bathy (ji,jj)=0.0 ; mbathy (ji,jj) = 0
1421	END IF
1422	END IF
1423	END DO
1424	END DO
1425	! ensure halo correct
1426	IF( lk_mpp ) THEN
1427	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1428	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1429	CALL lbc_lnk( risfdep, 'T', 1. )
1430	CALL lbc_lnk( bathy , 'T', 1. )
1431	ENDIF
1432	!
1433	! 2 Compute misfdep for ocean points (i.e. first wet level)
1434	! find the first ocean level such that the first level thickness
1435	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1436	! e3t_0 is the reference level thickness
1437	!
1438	! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf
1439	WHERE( risfdep(:,:) == 0.0 ) ; misfdep(:,:) = 1 ! open water or grounded : set misfdep to 1
1440	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
1441	END WHERE
1442	!
1443	DO jk = 2, jpkm1
1444	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1445	WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
1446	END DO
1447	!
1448	! 2.1 be sure misfdep not below mbathy
1449	! warning because of condition 4 we could have 'wet cell with misfdep below mbathy
1450	! risfdep of these cells will be fix later on (see 5)
1451	WHERE( misfdep > mbathy .AND. mbathy > 0 ) misfdep=mbathy
1452	!
1453	! 3.0 Assure 2 wet cells in the column at T point and along the edge.
1454	! first do T, then loop 4 times on U-1, U, V-1, V and then at T point
1455	! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
1456	!
1457	! find the deepest isfdep level that fit the 2 wet cell on the water column
1458	! on all the sides (still need 4 pass)
1459	DO jn = 1, 4
1460	zmisfdep = misfdep
1461	DO jj = 2, jpjm1
1462	DO ji = 2, jpim1
1463	! ISF cell only
1464	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1465	zmbathyij = mbathy(ji ,jj)
1466	! set ground edge value to jpk to skip it later on
1467	zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge
1468	zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge
1469	zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge
1470	zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge
1471	! shallowest bathy level
1472	zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1)
1473	! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the
1474	! water column
1475	jk = MIN(misfdep(ji,jj),zmbathyij-1)
1476	! update isf draft if needed
1477	IF ( jk < misfdep(ji,jj) ) THEN
1478	zmisfdep(ji,jj) = jk
1479	risfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1480	ELSE
1481	! isf depth not modify, nothing to do
1482	END IF
1483	ENDIF
1484	END DO
1485	END DO
1486	misfdep=zmisfdep
1487	! ensure halo correct before new pass
1488	IF( lk_mpp ) THEN
1489	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1490	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1491	CALL lbc_lnk( risfdep, 'T', 1. )
1492	CALL lbc_lnk( bathy , 'T', 1. )
1493	ENDIF
1494	END DO ! jn
1495	!
1496	! 3.1 condition block to assure connectivity every where beneath an ice shelf
1497	IF (ln_isfconnect) THEN
1498	zmask(:,:)=1.0
1499	WHERE (mbathy==0)
1500	zmask(:,:)=jpk
1501	END WHERE
1502
1503	zmbathy = mbathy
1504	zmisfdep = misfdep
1505	zrisfdep = risfdep
1506	WHERE (mbathy == 0) zmbathy=jpk
1507	DO jj = 2, jpjm1
1508	DO ji = 2, jpim1
1509	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1510	! what it should be
1511	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = should be connected
1512	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = should be connected
1513	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = should be connected
1514	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = should be connected
1515	! what it is
1516	umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk
1517	IF (misfdep(ji,jj) > zmbathy(ji+1,jj )) umskip1_r=1.0 ! 1 = no effective connection
1518	IF (misfdep(ji,jj) > zmbathy(ji-1,jj )) umskim1_r=1.0
1519	IF (misfdep(ji,jj) > zmbathy(ji ,jj+1)) vmskjp1_r=1.0
1520	IF (misfdep(ji,jj) > zmbathy(ji ,jj-1)) vmskjm1_r=1.0
1521	! defining level need for connectivity
1522	jk=NINT(MIN(zmbathy(ji+1,jj ) * umskip1_r * umskip1, &
1523	& zmbathy(ji-1,jj ) * umskim1_r * umskim1, &
1524	& zmbathy(ji ,jj+1) * vmskjp1_r * vmskjp1, &
1525	& zmbathy(ji ,jj-1) * vmskjm1_r * vmskjm1, &
1526	& jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on
1527	! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
1528	zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
1529	!
1530	! update isf draft if needed (need to be done here because next test is in meter and level)
1531	IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN
1532	jk = zmisfdep(ji,jj)
1533	zrisfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1534	END IF
1535
1536	! sanity check
1537	! check if we dig more than nn_kisfmax level or reach the surface
1538	! check if we dig more than rn_zisfmax meter
1539	! => if this is the case, undo what has been done before
1540	IF ( (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
1541	& .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj) )) ) THEN
1542	zmisfdep(ji,jj)=misfdep(ji,jj)
1543	zrisfdep(ji,jj)=risfdep(ji,jj)
1544	END IF
1545	END IF
1546	END DO
1547	END DO
1548	misfdep=zmisfdep
1549	risfdep=zrisfdep
1550	! ensure halo correct
1551	IF( lk_mpp ) THEN
1552	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1553	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1554	CALL lbc_lnk( risfdep, 'T', 1. )
1555	CALL lbc_lnk( bathy , 'T', 1. )
1556	ENDIF
1557	END IF
1558	!
1559	IF (ln_isfcheminey) THEN
1560	! 3.2 fill hole in ice shelf (ie cell with no velocity point)
1561	! => misfdep = MIN(misfdep at U, U-1, V, V-1)
1562	! risfdep = gdepw(misfdep) (ie full cell)
1563	zmisfdep = misfdep
1564	WHERE (mbathy == 0) zmisfdep=jpk ! grounded
1565	DO jj = 2, jpjm1
1566	DO ji = 2, jpim1
1567	ibtest = zmisfdep(ji ,jj )
1568	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
1569	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
1570	! case unconected wet cell (T point) where isfdraft(ii,jj) deeper
1571	! than bathy U-1/U+1/V-1/V+1
1572	IF( zmisfdep(ji,jj) > mbathy(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1573	IF( zmisfdep(ji,jj) > mbathy(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1574	IF( zmisfdep(ji,jj) > mbathy(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1575	IF( zmisfdep(ji,jj) > mbathy(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1576	! case unconected wet cell (T point) where bathy(ii,jj) shallower
1577	! than isfdraft U-1/U+1/V-1/V+1
1578	IF( mbathy(ji,jj) < zmisfdep(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1579	IF( mbathy(ji,jj) < zmisfdep(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1580	IF( mbathy(ji,jj) < zmisfdep(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1581	IF( mbathy(ji,jj) < zmisfdep(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1582	! if surround in fully mask => mask cell (1)
1583	! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)
1584	! if misfdep is not changed, nothing is done (the MAX is doing the
1585	! job)
1586	ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
1587	IF (misfdep(ji,jj) < ibtest) THEN
1588	misfdep(ji,jj) = ibtest
1589	risfdep(ji,jj) = gdepw_1d(ibtest)
1590	END IF
1591	ENDDO
1592	ENDDO
1593	WHERE( misfdep == jpk) ! ground case (1)
1594	mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp
1595	END WHERE
1596	!
1597	! ensure halo correct
1598	IF( lk_mpp ) THEN
1599	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1600	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1601	CALL lbc_lnk( risfdep, 'T', 1. )
1602	CALL lbc_lnk( bathy , 'T', 1. )
1603	END IF
1604	END IF
1605	!
1606	CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep)
1607	CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy)
1608	!
1609	IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')
1610
1611	END SUBROUTINE
1612
1613	SUBROUTINE zgr_sco
1614	!!----------------------------------------------------------------------
1615	!! * ROUTINE zgr_sco *
1616	!!
1617	!! ** Purpose : define the s-coordinate system
1618	!!
1619	!! ** Method : s-coordinate
1620	!! The depth of model levels is defined as the product of an
1621	!! analytical function by the local bathymetry, while the vertical
1622	!! scale factors are defined as the product of the first derivative
1623	!! of the analytical function by the bathymetry.
1624	!! (this solution save memory as depth and scale factors are not
1625	!! 3d fields)
1626	!! - Read bathymetry (in meters) at t-point and compute the
1627	!! bathymetry at u-, v-, and f-points.
1628	!! hbatu = mi( hbatt )
1629	!! hbatv = mj( hbatt )
1630	!! hbatf = mi( mj( hbatt ) )
1631	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1632	!! function and its derivative given as function.
1633	!! z_gsigt(k) = fssig (k )
1634	!! z_gsigw(k) = fssig (k-0.5)
1635	!! z_esigt(k) = fsdsig(k )
1636	!! z_esigw(k) = fsdsig(k-0.5)
1637	!! Three options for stretching are give, and they can be modified
1638	!! following the users requirements. Nevertheless, the output as
1639	!! well as the way to compute the model levels and scale factors
1640	!! must be respected in order to insure second order accuracy
1641	!! schemes.
1642	!!
1643	!! The three methods for stretching available are:
1644	!!
1645	!! s_sh94 (Song and Haidvogel 1994)
1646	!! a sinh/tanh function that allows sigma and stretched sigma
1647	!!
1648	!! s_sf12 (Siddorn and Furner 2012?)
1649	!! allows the maintenance of fixed surface and or
1650	!! bottom cell resolutions (cf. geopotential coordinates)
1651	!! within an analytically derived stretched S-coordinate framework.
1652	!!
1653	!! s_tanh (Madec et al 1996)
1654	!! a cosh/tanh function that gives stretched coordinates
1655	!!
1656	!!----------------------------------------------------------------------
1657	!
1658	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1659	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1660	INTEGER :: ios ! Local integer output status for namelist read
1661	REAL(wp) :: zrmax, ztaper ! temporary scalars
1662	REAL(wp) :: zrfact
1663	!
1664	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1665	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1666
1667	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1668	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1669	!!----------------------------------------------------------------------
1670	!
1671	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1672	!
1673	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1674	!
1675	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1676	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1677	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1678
1679	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1680	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1681	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1682	IF(lwm) WRITE ( numond, namzgr_sco )
1683
1684	IF(lwp) THEN ! control print
1685	WRITE(numout,*)
1686	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1687	WRITE(numout,*) '~~~~~~~~~~~'
1688	WRITE(numout,*) ' Namelist namzgr_sco'
1689	WRITE(numout,*) ' stretching coeffs '
1690	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1691	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1692	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1693	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1694	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1695	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1696	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1697	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1698	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1699	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1700	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1701	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1702	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1703	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1704	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1705	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1706	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1707	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1708	ENDIF
1709
1710	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1711	hifu(:,:) = rn_sbot_min
1712	hifv(:,:) = rn_sbot_min
1713	hiff(:,:) = rn_sbot_min
1714
1715	! ! set maximum ocean depth
1716	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1717
1718	DO jj = 1, jpj
1719	DO ji = 1, jpi
1720	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1721	END DO
1722	END DO
1723	! ! =============================
1724	! ! Define the envelop bathymetry (hbatt)
1725	! ! =============================
1726	! use r-value to create hybrid coordinates
1727	zenv(:,:) = bathy(:,:)
1728	!
1729	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1730	DO jj = 1, jpj
1731	DO ji = 1, jpi
1732	IF( bathy(ji,jj) == 0._wp ) THEN
1733	iip1 = MIN( ji+1, jpi )
1734	ijp1 = MIN( jj+1, jpj )
1735	iim1 = MAX( ji-1, 1 )
1736	ijm1 = MAX( jj-1, 1 )
1737	IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
1738	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
1739	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) THEN
1740	zenv(ji,jj) = rn_sbot_min
1741	ENDIF
1742	ENDIF
1743	END DO
1744	END DO
1745	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1746	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1747	!
1748	! smooth the bathymetry (if required)
1749	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1750	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1751	!
1752	jl = 0
1753	zrmax = 1._wp
1754	!
1755	!
1756	! set scaling factor used in reducing vertical gradients
1757	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1758	!
1759	! initialise temporary evelope depth arrays
1760	ztmpi1(:,:) = zenv(:,:)
1761	ztmpi2(:,:) = zenv(:,:)
1762	ztmpj1(:,:) = zenv(:,:)
1763	ztmpj2(:,:) = zenv(:,:)
1764	!
1765	! initialise temporary r-value arrays
1766	zri(:,:) = 1._wp
1767	zrj(:,:) = 1._wp
1768	! ! ================ !
1769	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1770	! ! ================ !
1771	jl = jl + 1
1772	zrmax = 0._wp
1773	! we set zrmax from previous r-values (zri and zrj) first
1774	! if set after current r-value calculation (as previously)
1775	! we could exit DO WHILE prematurely before checking r-value
1776	! of current zenv
1777	DO jj = 1, nlcj
1778	DO ji = 1, nlci
1779	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1780	END DO
1781	END DO
1782	zri(:,:) = 0._wp
1783	zrj(:,:) = 0._wp
1784	DO jj = 1, nlcj
1785	DO ji = 1, nlci
1786	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1787	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1788	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1789	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1790	END IF
1791	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1792	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1793	END IF
1794	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
1795	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
1796	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
1797	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
1798	END DO
1799	END DO
1800	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1801	!
1802	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
1803	!
1804	DO jj = 1, nlcj
1805	DO ji = 1, nlci
1806	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1807	END DO
1808	END DO
1809	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1810	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1811	! ! ================ !
1812	END DO ! End loop !
1813	! ! ================ !
1814	DO jj = 1, jpj
1815	DO ji = 1, jpi
1816	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1817	END DO
1818	END DO
1819	!
1820	! Envelope bathymetry saved in hbatt
1821	hbatt(:,:) = zenv(:,:)
1822	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1823	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1824	DO jj = 1, jpj
1825	DO ji = 1, jpi
1826	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1827	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1828	END DO
1829	END DO
1830	ENDIF
1831	!
1832	IF(lwp) THEN ! Control print
1833	WRITE(numout,*)
1834	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1835	WRITE(numout,*)
1836	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1837	IF( nprint == 1 ) THEN
1838	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1839	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1840	ENDIF
1841	ENDIF
1842
1843	! ! ==============================
1844	! ! hbatu, hbatv, hbatf fields
1845	! ! ==============================
1846	IF(lwp) THEN
1847	WRITE(numout,*)
1848	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1849	ENDIF
1850	hbatu(:,:) = rn_sbot_min
1851	hbatv(:,:) = rn_sbot_min
1852	hbatf(:,:) = rn_sbot_min
1853	DO jj = 1, jpjm1
1854	DO ji = 1, jpim1 ! NO vector opt.
1855	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1856	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1857	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1858	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1859	END DO
1860	END DO
1861	!
1862	! Apply lateral boundary condition
1863	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1864	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
1865	DO jj = 1, jpj
1866	DO ji = 1, jpi
1867	IF( hbatu(ji,jj) == 0._wp ) THEN
1868	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1869	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1870	ENDIF
1871	END DO
1872	END DO
1873	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
1874	DO jj = 1, jpj
1875	DO ji = 1, jpi
1876	IF( hbatv(ji,jj) == 0._wp ) THEN
1877	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1878	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1879	ENDIF
1880	END DO
1881	END DO
1882	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
1883	DO jj = 1, jpj
1884	DO ji = 1, jpi
1885	IF( hbatf(ji,jj) == 0._wp ) THEN
1886	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1887	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1888	ENDIF
1889	END DO
1890	END DO
1891
1892	!!bug: key_helsinki a verifer
1893	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1894	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1895	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1896	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1897
1898	IF( nprint == 1 .AND. lwp ) THEN
1899	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1900	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1901	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1902	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1903	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1904	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1905	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1906	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1907	ENDIF
1908	!! helsinki
1909
1910	! ! =======================
1911	! ! s-ccordinate fields (gdep., e3.)
1912	! ! =======================
1913	!
1914	! non-dimensional "sigma" for model level depth at w- and t-levels
1915
1916
1917	!========================================================================
1918	! Song and Haidvogel 1994 (ln_s_sh94=T)
1919	! Siddorn and Furner 2012 (ln_sf12=T)
1920	! or tanh function (both false)
1921	!========================================================================
1922	IF ( ln_s_sh94 ) THEN
1923	CALL s_sh94()
1924	ELSE IF ( ln_s_sf12 ) THEN
1925	CALL s_sf12()
1926	ELSE
1927	CALL s_tanh()
1928	ENDIF
1929
1930	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
1931	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
1932	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
1933	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
1934	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
1935	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
1936	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1937
1938	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1939	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1940	!
1941	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
1942	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
1943	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
1944	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
1945	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
1946	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
1947	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
1948
1949	#if defined key_agrif
1950	! Ensure meaningful vertical scale factors in ghost lines/columns
1951	IF( .NOT. Agrif_Root() ) THEN
1952	!
1953	IF((nbondi == -1).OR.(nbondi == 2)) THEN
1954	e3u_0(1,:,:) = e3u_0(2,:,:)
1955	ENDIF
1956	!
1957	IF((nbondi == 1).OR.(nbondi == 2)) THEN
1958	e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
1959	ENDIF
1960	!
1961	IF((nbondj == -1).OR.(nbondj == 2)) THEN
1962	e3v_0(:,1,:) = e3v_0(:,2,:)
1963	ENDIF
1964	!
1965	IF((nbondj == 1).OR.(nbondj == 2)) THEN
1966	e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
1967	ENDIF
1968	!
1969	ENDIF
1970	#endif
1971
1972	fsdept(:,:,:) = gdept_0 (:,:,:)
1973	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1974	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1975	fse3t (:,:,:) = e3t_0 (:,:,:)
1976	fse3u (:,:,:) = e3u_0 (:,:,:)
1977	fse3v (:,:,:) = e3v_0 (:,:,:)
1978	fse3f (:,:,:) = e3f_0 (:,:,:)
1979	fse3w (:,:,:) = e3w_0 (:,:,:)
1980	fse3uw(:,:,:) = e3uw_0 (:,:,:)
1981	fse3vw(:,:,:) = e3vw_0 (:,:,:)
1982	!!
1983	! HYBRID :
1984	DO jj = 1, jpj
1985	DO ji = 1, jpi
1986	DO jk = 1, jpkm1
1987	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1988	END DO
1989	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
1990	END DO
1991	END DO
1992	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1993	& ' MAX ', MAXVAL( mbathy(:,:) )
1994
1995	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1996	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1997	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
1998	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
1999	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
2000	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
2001	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
2002	& ' w ', MINVAL( e3w_0 (:,:,:) )
2003
2004	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
2005	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
2006	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
2007	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
2008	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
2009	& ' w ', MAXVAL( e3w_0 (:,:,:) )
2010	ENDIF
2011	! END DO
2012	IF(lwp) THEN ! selected vertical profiles
2013	WRITE(numout,*)
2014	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2015	WRITE(numout,*) ' ~~~~~~ --------------------'
2016	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2017	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
2018	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2019	DO jj = mj0(20), mj1(20)
2020	DO ji = mi0(20), mi1(20)
2021	WRITE(numout,*)
2022	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2023	WRITE(numout,*) ' ~~~~~~ --------------------'
2024	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2025	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2026	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2027	END DO
2028	END DO
2029	DO jj = mj0(74), mj1(74)
2030	DO ji = mi0(100), mi1(100)
2031	WRITE(numout,*)
2032	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2033	WRITE(numout,*) ' ~~~~~~ --------------------'
2034	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2035	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2036	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2037	END DO
2038	END DO
2039	ENDIF
2040
2041	!================================================================================
2042	! check the coordinate makes sense
2043	!================================================================================
2044	DO ji = 1, jpi
2045	DO jj = 1, jpj
2046
2047	IF( hbatt(ji,jj) > 0._wp) THEN
2048	DO jk = 1, mbathy(ji,jj)
2049	! check coordinate is monotonically increasing
2050	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2051	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2052	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2053	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2054	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2055	CALL ctl_stop( ctmp1 )
2056	ENDIF
2057	! and check it has never gone negative
2058	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2059	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2060	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2061	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2062	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2063	CALL ctl_stop( ctmp1 )
2064	ENDIF
2065	! and check it never exceeds the total depth
2066	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2067	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2068	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2069	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2070	CALL ctl_stop( ctmp1 )
2071	ENDIF
2072	END DO
2073
2074	DO jk = 1, mbathy(ji,jj)-1
2075	! and check it never exceeds the total depth
2076	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2077	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2078	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2079	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2080	CALL ctl_stop( ctmp1 )
2081	ENDIF
2082	END DO
2083
2084	ENDIF
2085
2086	END DO
2087	END DO
2088	!
2089	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2090	!
2091	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
2092	!
2093	END SUBROUTINE zgr_sco
2094
2095	!!======================================================================
2096	SUBROUTINE s_sh94()
2097
2098	!!----------------------------------------------------------------------
2099	!! * ROUTINE s_sh94 *
2100	!!
2101	!! ** Purpose : stretch the s-coordinate system
2102	!!
2103	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
2104	!! mixed S/sigma coordinate
2105	!!
2106	!! Reference : Song and Haidvogel 1994.
2107	!!----------------------------------------------------------------------
2108	!
2109	INTEGER :: ji, jj, jk ! dummy loop argument
2110	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2111	!
2112	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2113	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2114
2115	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2116	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2117
2118	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2119	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2120	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2121	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2122
2123	DO ji = 1, jpi
2124	DO jj = 1, jpj
2125
2126	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2127	DO jk = 1, jpk
2128	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2129	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2130	END DO
2131	ELSE ! shallow water, uniform sigma
2132	DO jk = 1, jpk
2133	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
2134	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2135	END DO
2136	ENDIF
2137	!
2138	DO jk = 1, jpkm1
2139	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2140	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2141	END DO
2142	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
2143	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2144	!
2145	! Coefficients for vertical depth as the sum of e3w scale factors
2146	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2147	DO jk = 2, jpk
2148	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2149	END DO
2150	!
2151	DO jk = 1, jpk
2152	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2153	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2154	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
2155	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
2156	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
2157	END DO
2158	!
2159	END DO ! for all jj's
2160	END DO ! for all ji's
2161
2162	DO ji = 1, jpim1
2163	DO jj = 1, jpjm1
2164	DO jk = 1, jpk
2165	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
2166	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2167	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
2168	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2169	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
2170	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
2171	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2172	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
2173	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2174	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
2175	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2176	!
2177	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2178	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2179	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2180	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2181	!
2182	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2183	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2184	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2185	END DO
2186	END DO
2187	END DO
2188
2189	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2190	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2191
2192	END SUBROUTINE s_sh94
2193
2194	SUBROUTINE s_sf12
2195
2196	!!----------------------------------------------------------------------
2197	!! * ROUTINE s_sf12 *
2198	!!
2199	!! ** Purpose : stretch the s-coordinate system
2200	!!
2201	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
2202	!! mixed S/sigma/Z coordinate
2203	!!
2204	!! This method allows the maintenance of fixed surface and or
2205	!! bottom cell resolutions (cf. geopotential coordinates)
2206	!! within an analytically derived stretched S-coordinate framework.
2207	!!
2208	!!
2209	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2210	!!----------------------------------------------------------------------
2211	!
2212	INTEGER :: ji, jj, jk ! dummy loop argument
2213	REAL(wp) :: zsmth ! smoothing around critical depth
2214	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
2215	!
2216	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2217	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2218
2219	!
2220	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2221	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2222
2223	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2224	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2225	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2226	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2227
2228	DO ji = 1, jpi
2229	DO jj = 1, jpj
2230
2231	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2232
2233	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
2234	! could be changed by users but care must be taken to do so carefully
2235	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2236
2237	zzs = rn_zs / hbatt(ji,jj)
2238
2239	IF (rn_efold /= 0.0_wp) THEN
2240	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2241	ELSE
2242	zsmth = 1.0_wp
2243	ENDIF
2244
2245	DO jk = 1, jpk
2246	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2247	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2248	ENDDO
2249	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
2250	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
2251
2252	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2253
2254	DO jk = 1, jpk
2255	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2256	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2257	END DO
2258
2259	ELSE ! shallow water, z coordinates
2260
2261	DO jk = 1, jpk
2262	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2263	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2264	END DO
2265
2266	ENDIF
2267
2268	DO jk = 1, jpkm1
2269	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2270	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2271	END DO
2272	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
2273	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2274
2275	! Coefficients for vertical depth as the sum of e3w scale factors
2276	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2277	DO jk = 2, jpk
2278	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2279	END DO
2280
2281	DO jk = 1, jpk
2282	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2283	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2284	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2285	END DO
2286
2287	ENDDO ! for all jj's
2288	ENDDO ! for all ji's
2289
2290	DO ji=1,jpi-1
2291	DO jj=1,jpj-1
2292
2293	DO jk = 1, jpk
2294	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
2295	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2296	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
2297	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2298	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
2299	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
2300	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2301	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
2302	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2303	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
2304	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2305
2306	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2307	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2308	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2309	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2310	!
2311	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2312	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2313	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2314	END DO
2315
2316	ENDDO
2317	ENDDO
2318	!
2319	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2320	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2321	CALL lbc_lnk(e3w_0 ,'T',1.)
2322	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2323	!
2324	! ! =============
2325
2326	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2327	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2328
2329	END SUBROUTINE s_sf12
2330
2331	SUBROUTINE s_tanh()
2332
2333	!!----------------------------------------------------------------------
2334	!! * ROUTINE s_tanh*
2335	!!
2336	!! ** Purpose : stretch the s-coordinate system
2337	!!
2338	!! ** Method : s-coordinate stretch
2339	!!
2340	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2341	!!----------------------------------------------------------------------
2342
2343	INTEGER :: ji, jj, jk ! dummy loop argument
2344	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2345
2346	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2347	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2348
2349	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2350	CALL wrk_alloc( jpk, z_esigt, z_esigw )
2351
2352	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
2353	z_esigt = 0._wp ; z_esigw = 0._wp
2354
2355	DO jk = 1, jpk
2356	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2357	z_gsigt(jk) = -fssig( REAL(jk,wp) )
2358	END DO
2359	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
2360	!
2361	! Coefficients for vertical scale factors at w-, t- levels
2362	!!gm bug : define it from analytical function, not like juste bellow....
2363	!!gm or betteroffer the 2 possibilities....
2364	DO jk = 1, jpkm1
2365	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
2366	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2367	END DO
2368	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
2369	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2370	!
2371	! Coefficients for vertical depth as the sum of e3w scale factors
2372	z_gsi3w(1) = 0.5_wp * z_esigw(1)
2373	DO jk = 2, jpk
2374	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2375	END DO
2376	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2377	DO jk = 1, jpk
2378	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2379	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2380	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
2381	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
2382	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
2383	END DO
2384	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2385	DO jj = 1, jpj
2386	DO ji = 1, jpi
2387	DO jk = 1, jpk
2388	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2389	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2390	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2391	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2392	!
2393	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2394	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2395	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2396	END DO
2397	END DO
2398	END DO
2399
2400	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2401	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2402
2403	END SUBROUTINE s_tanh
2404
2405	FUNCTION fssig( pk ) RESULT( pf )
2406	!!----------------------------------------------------------------------
2407	!! * ROUTINE fssig *
2408	!!
2409	!! ** Purpose : provide the analytical function in s-coordinate
2410	!!
2411	!! ** Method : the function provide the non-dimensional position of
2412	!! T and W (i.e. between 0 and 1)
2413	!! T-points at integer values (between 1 and jpk)
2414	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2415	!!----------------------------------------------------------------------
2416	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2417	REAL(wp) :: pf ! sigma value
2418	!!----------------------------------------------------------------------
2419	!
2420	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2421	& - TANH( rn_thetb * rn_theta ) ) &
2422	& * ( COSH( rn_theta ) &
2423	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2424	& / ( 2._wp * SINH( rn_theta ) )
2425	!
2426	END FUNCTION fssig
2427
2428
2429	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2430	!!----------------------------------------------------------------------
2431	!! * ROUTINE fssig1 *
2432	!!
2433	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2434	!!
2435	!! ** Method : the function provides the non-dimensional position of
2436	!! T and W (i.e. between 0 and 1)
2437	!! T-points at integer values (between 1 and jpk)
2438	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2439	!!----------------------------------------------------------------------
2440	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2441	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2442	REAL(wp) :: pf1 ! sigma value
2443	!!----------------------------------------------------------------------
2444	!
2445	IF ( rn_theta == 0 ) then ! uniform sigma
2446	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2447	ELSE ! stretched sigma
2448	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2449	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2450	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2451	ENDIF
2452	!
2453	END FUNCTION fssig1
2454
2455
2456	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2457	!!----------------------------------------------------------------------
2458	!! * ROUTINE fgamma *
2459	!!
2460	!! ** Purpose : provide analytical function for the s-coordinate
2461	!!
2462	!! ** Method : the function provides the non-dimensional position of
2463	!! T and W (i.e. between 0 and 1)
2464	!! T-points at integer values (between 1 and jpk)
2465	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2466	!!
2467	!! This method allows the maintenance of fixed surface and or
2468	!! bottom cell resolutions (cf. geopotential coordinates)
2469	!! within an analytically derived stretched S-coordinate framework.
2470	!!
2471	!! Reference : Siddorn and Furner, in prep
2472	!!----------------------------------------------------------------------
2473	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2474	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2475	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2476	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2477	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2478	REAL(wp) :: za1,za2,za3 ! local variables
2479	REAL(wp) :: zn1,zn2 ! local variables
2480	REAL(wp) :: za,zb,zx ! local variables
2481	integer :: jk
2482	!!----------------------------------------------------------------------
2483	!
2484
2485	zn1 = 1./(jpk-1.)
2486	zn2 = 1. - zn1
2487
2488	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2489	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2490	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2491
2492	za = pzb - za3*(pzs-za1)-za2
2493	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2494	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2495	zx = 1.0_wp-za/2.0_wp-zb
2496
2497	DO jk = 1, jpk
2498	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2499	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2500	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2501	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2502	ENDDO
2503
2504	!
2505	END FUNCTION fgamma
2506
2507	!!======================================================================
2508	END MODULE domzgr

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

Download in other formats: