Context Navigation

domzgr.F90 @ 11248

Last change on this file since 11248 was 11248, checked in by mathiot, 5 years ago
add condition to fill isolated grid point in the bathymetry in zgr_isf to ensure compatibility with ice shelf draft and properly mask bathy,risfdep,mbathy,misfdep after filling subglacial lakes
File size: 126.7 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
7	!! ! 1997-07 (G. Madec) lbc_lnk call
8	!! ! 1997-04 (J.-O. Beismann)
9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
19	!! 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 :: ios, icompt
1290	INTEGER :: ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
1291	INTEGER :: zmbathyij, zmbathyip1, zmbathyim1, zmbathyjp1, zmbathyjm1
1292	INTEGER , POINTER, DIMENSION(:,:) :: zmisfdep, zmbathy ! 2D workspace (ISH)
1293	!
1294	REAL(wp) :: zdepth, vmskjp1, vmskjm1, umskip1, umskim1, umskip1_r, umskim1_r, vmskjp1_r, vmskjm1_r
1295	REAL(wp) :: zisfdep_min
1296	REAL(wp), POINTER, DIMENSION(:,:) :: zdummy, zmask, zrisfdep ! 2D workspace (ISH)
1297	!
1298	!!-----------------
1299	! NAMELIST
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	!
1426	! ensure halo correct
1427	IF( lk_mpp ) THEN
1428	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1429	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1430	CALL lbc_lnk( risfdep, 'T', 1. )
1431	CALL lbc_lnk( bathy , 'T', 1. )
1432	ENDIF
1433	!
1434	! 2 Compute misfdep for ocean points (i.e. first wet level)
1435	! find the first ocean level such that the first level thickness
1436	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1437	! e3t_0 is the reference level thickness
1438	!
1439	! 2.0 set misfdep to 1 on open water and initialisation beneath ice shelf
1440	WHERE( risfdep(:,:) == 0.0 ) ; misfdep(:,:) = 1 ! open water or grounded : set misfdep to 1
1441	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
1442	END WHERE
1443	!
1444	DO jk = 2, jpkm1
1445	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1446	WHERE( risfdep(:,:) > 0.0 .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
1447	END DO
1448	!
1449	! 2.1 fill isolated grid point in the bathymetry
1450	! need to be done here to fix compatibility (will be done again in zgr_bat_ctl)
1451	DO jj = 2, jpjm1
1452	DO ji = 2, jpim1
1453	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
1454	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
1455	IF( ibtest < mbathy(ji,jj) ) THEN
1456	mbathy(ji,jj) = ibtest
1457	icompt = icompt + 1
1458	END IF
1459	END DO
1460	END DO
1461	!
1462	! ensure halo correct
1463	zdummy(:,:) = FLOAT( mbathy(:,:) ) ; CALL lbc_lnk( zdummy, 'T', 1._wp ) ; mbathy(:,:) = INT( zdummy(:,:) )
1464	!
1465	IF( lk_mpp ) CALL mpp_sum( icompt )
1466	IF( icompt == 0 ) THEN
1467	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
1468	ELSE
1469	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
1470	ENDIF
1471	!
1472	! 2.2 be sure misfdep not below mbathy
1473	! warning because of condition 2.0 and 2.1 we could have 'wet cell with misfdep below mbathy
1474	! risfdep of these cells will be fix later on
1475	WHERE( misfdep > mbathy ) misfdep(:,:) = MAX( 1, mbathy(:,:) )
1476	!
1477	! 3.0 Assure 2 wet cells in the column at T point and along the edge.
1478	! first do T, then loop 4 times on U-1, U, V-1, V and then at T point
1479	! If un-luky, digging cell U-1 can trigger case for U+1, then V-1, then V+1
1480	!
1481	! find the deepest isfdep level that fit the 2 wet cell on the water column
1482	! on all the sides (still need 4 pass)
1483	DO jn = 1, 4
1484	zmisfdep = misfdep
1485	DO jj = 2, jpjm1
1486	DO ji = 2, jpim1
1487	! ISF cell only
1488	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1489	zmbathyij = mbathy(ji ,jj)
1490	! set ground edge value to jpk to skip it later on
1491	zmbathyip1 = mbathy(ji+1,jj) ; IF (zmbathyip1 < misfdep(ji,jj)) zmbathyip1 = jpk ! not wet cell in common on this edge
1492	zmbathyim1 = mbathy(ji-1,jj) ; IF (zmbathyim1 < misfdep(ji,jj)) zmbathyim1 = jpk ! not wet cell in common on this edge
1493	zmbathyjp1 = mbathy(ji,jj+1) ; IF (zmbathyjp1 < misfdep(ji,jj)) zmbathyjp1 = jpk ! not wet cell in common on this edge
1494	zmbathyjm1 = mbathy(ji,jj-1) ; IF (zmbathyjm1 < misfdep(ji,jj)) zmbathyjm1 = jpk ! not wet cell in common on this edge
1495	! shallowest bathy level
1496	zmbathyij = MIN(zmbathyij,zmbathyip1,zmbathyim1,zmbathyjp1,zmbathyjm1)
1497	! misfdep need to be <= zmbathyij-1 to fit 2 wet cell on the
1498	! water column
1499	jk = MIN(misfdep(ji,jj),zmbathyij-1)
1500	! update isf draft if needed
1501	IF ( jk < misfdep(ji,jj) ) THEN
1502	zmisfdep(ji,jj) = jk
1503	risfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1504	ELSE
1505	! isf depth not modify, nothing to do
1506	END IF
1507	ENDIF
1508	END DO
1509	END DO
1510	misfdep=zmisfdep
1511	! ensure halo correct before new pass
1512	IF( lk_mpp ) THEN
1513	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1514	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1515	CALL lbc_lnk( risfdep, 'T', 1. )
1516	CALL lbc_lnk( bathy , 'T', 1. )
1517	ENDIF
1518	END DO ! jn
1519	!
1520	! 3.1 condition block to assure connectivity every where beneath an ice shelf
1521	IF (ln_isfconnect) THEN
1522	zmask(:,:)=1.0
1523	WHERE (mbathy==0)
1524	zmask(:,:)=jpk
1525	END WHERE
1526
1527	zmbathy = mbathy
1528	zmisfdep = misfdep
1529	zrisfdep = risfdep
1530	WHERE (mbathy == 0) zmbathy=jpk
1531	DO jj = 2, jpjm1
1532	DO ji = 2, jpim1
1533	IF( (misfdep(ji,jj) > 1) .AND. (mbathy(ji,jj) > 0) ) THEN
1534	! what it should be
1535	umskip1=zmask(ji,jj)*zmask(ji+1,jj ) ! 1 = should be connected
1536	umskim1=zmask(ji,jj)*zmask(ji-1,jj ) ! 1 = should be connected
1537	vmskjp1=zmask(ji,jj)*zmask(ji ,jj+1) ! 1 = should be connected
1538	vmskjm1=zmask(ji,jj)*zmask(ji ,jj-1) ! 1 = should be connected
1539	! what it is
1540	umskip1_r=jpk ; umskim1_r=jpk; vmskjp1_r=jpk; vmskjm1_r=jpk
1541	IF (misfdep(ji,jj) > zmbathy(ji+1,jj )) umskip1_r=1.0 ! 1 = no effective connection
1542	IF (misfdep(ji,jj) > zmbathy(ji-1,jj )) umskim1_r=1.0
1543	IF (misfdep(ji,jj) > zmbathy(ji ,jj+1)) vmskjp1_r=1.0
1544	IF (misfdep(ji,jj) > zmbathy(ji ,jj-1)) vmskjm1_r=1.0
1545	! defining level need for connectivity
1546	jk=NINT(MIN(zmbathy(ji+1,jj ) * umskip1_r * umskip1, &
1547	& zmbathy(ji-1,jj ) * umskim1_r * umskim1, &
1548	& zmbathy(ji ,jj+1) * vmskjp1_r * vmskjp1, &
1549	& zmbathy(ji ,jj-1) * vmskjm1_r * vmskjm1, &
1550	& jpk+1.0 )) ! add jpk in the MIN to avoid out of boundary later on
1551	! if connectivity is OK or no connection needed (grounding line) or grounded, zmisfdep=misfdep
1552	zmisfdep(ji,jj)=MIN(misfdep(ji,jj),jk-1)
1553	!
1554	! update isf draft if needed (need to be done here because next test is in meter and level)
1555	IF ( zmisfdep(ji,jj) < misfdep(ji,jj) ) THEN
1556	jk = zmisfdep(ji,jj)
1557	zrisfdep(ji,jj) = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1558	END IF
1559
1560	! sanity check
1561	! check if we dig more than nn_kisfmax level or reach the surface
1562	! check if we dig more than rn_zisfmax meter
1563	! => if this is the case, undo what has been done before
1564	IF ( (misfdep(ji,jj)-zmisfdep(ji,jj) > MIN(nn_kisfmax,misfdep(ji,jj)-2)) &
1565	& .OR. (risfdep(ji,jj)-zrisfdep(ji,jj) > MIN(rn_zisfmax,risfdep(ji,jj) )) ) THEN
1566	zmisfdep(ji,jj)=misfdep(ji,jj)
1567	zrisfdep(ji,jj)=risfdep(ji,jj)
1568	END IF
1569	END IF
1570	END DO
1571	END DO
1572	misfdep=zmisfdep
1573	risfdep=zrisfdep
1574	! ensure halo correct
1575	IF( lk_mpp ) THEN
1576	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1577	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1578	CALL lbc_lnk( risfdep, 'T', 1. )
1579	CALL lbc_lnk( bathy , 'T', 1. )
1580	ENDIF
1581	END IF
1582	!
1583	IF (ln_isfcheminey) THEN
1584	! 3.2 fill hole in ice shelf (ie cell with no velocity point)
1585	! => misfdep = MIN(misfdep at U, U-1, V, V-1)
1586	! risfdep = gdepw(misfdep) (ie full cell)
1587	zmisfdep = misfdep
1588	WHERE (mbathy == 0) zmisfdep=jpk ! grounded
1589	DO jj = 2, jpjm1
1590	DO ji = 2, jpim1
1591	ibtest = zmisfdep(ji ,jj )
1592	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
1593	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
1594	! case unconected wet cell (T point) where isfdraft(ii,jj) deeper
1595	! than bathy U-1/U+1/V-1/V+1
1596	IF( zmisfdep(ji,jj) > mbathy(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1597	IF( zmisfdep(ji,jj) > mbathy(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1598	IF( zmisfdep(ji,jj) > mbathy(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1599	IF( zmisfdep(ji,jj) > mbathy(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1600	! case unconected wet cell (T point) where bathy(ii,jj) shallower
1601	! than isfdraft U-1/U+1/V-1/V+1
1602	IF( mbathy(ji,jj) < zmisfdep(ji-1,jj ) ) ibtestim1 = jpk ! mask at U-1
1603	IF( mbathy(ji,jj) < zmisfdep(ji+1,jj ) ) ibtestip1 = jpk ! mask at U+1
1604	IF( mbathy(ji,jj) < zmisfdep(ji ,jj-1) ) ibtestjm1 = jpk ! mask at V-1
1605	IF( mbathy(ji,jj) < zmisfdep(ji ,jj+1) ) ibtestjp1 = jpk ! mask at V+1
1606	! if surround in fully mask => mask cell (1)
1607	! if hole in the ice shelf, set to the min of surrounding (the MIN is doing the job)
1608	! if misfdep is not changed, nothing is done (the MAX is doing the
1609	! job)
1610	ibtest = MAX(ibtest,MIN(ibtestim1, ibtestip1, ibtestjm1,ibtestjp1))
1611	IF (misfdep(ji,jj) < ibtest) THEN
1612	misfdep(ji,jj) = ibtest
1613	risfdep(ji,jj) = gdepw_1d(ibtest)
1614	END IF
1615	ENDDO
1616	ENDDO
1617	WHERE( misfdep == jpk) ! ground case (1)
1618	mbathy = 0 ; bathy = 0.0_wp ; misfdep = 1 ; risfdep = 0.0_wp
1619	END WHERE
1620	!
1621	! ensure halo correct
1622	IF( lk_mpp ) THEN
1623	zdummy(:,:) = FLOAT( misfdep(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); misfdep(:,:) = INT( zdummy(:,:) )
1624	zdummy(:,:) = FLOAT( mbathy(:,:) ); CALL lbc_lnk( zdummy, 'T', 1. ); mbathy(:,:) = INT( zdummy(:,:) )
1625	CALL lbc_lnk( risfdep, 'T', 1. )
1626	CALL lbc_lnk( bathy , 'T', 1. )
1627	END IF
1628	END IF
1629	!
1630	CALL wrk_dealloc( jpi, jpj, zmask, zdummy, zrisfdep)
1631	CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy)
1632	!
1633	IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')
1634
1635	END SUBROUTINE
1636
1637	SUBROUTINE zgr_sco
1638	!!----------------------------------------------------------------------
1639	!! * ROUTINE zgr_sco *
1640	!!
1641	!! ** Purpose : define the s-coordinate system
1642	!!
1643	!! ** Method : s-coordinate
1644	!! The depth of model levels is defined as the product of an
1645	!! analytical function by the local bathymetry, while the vertical
1646	!! scale factors are defined as the product of the first derivative
1647	!! of the analytical function by the bathymetry.
1648	!! (this solution save memory as depth and scale factors are not
1649	!! 3d fields)
1650	!! - Read bathymetry (in meters) at t-point and compute the
1651	!! bathymetry at u-, v-, and f-points.
1652	!! hbatu = mi( hbatt )
1653	!! hbatv = mj( hbatt )
1654	!! hbatf = mi( mj( hbatt ) )
1655	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1656	!! function and its derivative given as function.
1657	!! z_gsigt(k) = fssig (k )
1658	!! z_gsigw(k) = fssig (k-0.5)
1659	!! z_esigt(k) = fsdsig(k )
1660	!! z_esigw(k) = fsdsig(k-0.5)
1661	!! Three options for stretching are give, and they can be modified
1662	!! following the users requirements. Nevertheless, the output as
1663	!! well as the way to compute the model levels and scale factors
1664	!! must be respected in order to insure second order accuracy
1665	!! schemes.
1666	!!
1667	!! The three methods for stretching available are:
1668	!!
1669	!! s_sh94 (Song and Haidvogel 1994)
1670	!! a sinh/tanh function that allows sigma and stretched sigma
1671	!!
1672	!! s_sf12 (Siddorn and Furner 2012?)
1673	!! allows the maintenance of fixed surface and or
1674	!! bottom cell resolutions (cf. geopotential coordinates)
1675	!! within an analytically derived stretched S-coordinate framework.
1676	!!
1677	!! s_tanh (Madec et al 1996)
1678	!! a cosh/tanh function that gives stretched coordinates
1679	!!
1680	!!----------------------------------------------------------------------
1681	!
1682	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1683	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1684	INTEGER :: ios ! Local integer output status for namelist read
1685	REAL(wp) :: zrmax, ztaper ! temporary scalars
1686	REAL(wp) :: zrfact
1687	!
1688	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1689	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1690
1691	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1692	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1693	!!----------------------------------------------------------------------
1694	!
1695	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1696	!
1697	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1698	!
1699	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1700	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1701	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1702
1703	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1704	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1705	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1706	IF(lwm) WRITE ( numond, namzgr_sco )
1707
1708	IF(lwp) THEN ! control print
1709	WRITE(numout,*)
1710	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1711	WRITE(numout,*) '~~~~~~~~~~~'
1712	WRITE(numout,*) ' Namelist namzgr_sco'
1713	WRITE(numout,*) ' stretching coeffs '
1714	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1715	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1716	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1717	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1718	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1719	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1720	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1721	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1722	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1723	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1724	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1725	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1726	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1727	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1728	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1729	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1730	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1731	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1732	ENDIF
1733
1734	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1735	hifu(:,:) = rn_sbot_min
1736	hifv(:,:) = rn_sbot_min
1737	hiff(:,:) = rn_sbot_min
1738
1739	! ! set maximum ocean depth
1740	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1741
1742	DO jj = 1, jpj
1743	DO ji = 1, jpi
1744	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1745	END DO
1746	END DO
1747	! ! =============================
1748	! ! Define the envelop bathymetry (hbatt)
1749	! ! =============================
1750	! use r-value to create hybrid coordinates
1751	zenv(:,:) = bathy(:,:)
1752	!
1753	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1754	DO jj = 1, jpj
1755	DO ji = 1, jpi
1756	IF( bathy(ji,jj) == 0._wp ) THEN
1757	iip1 = MIN( ji+1, jpi )
1758	ijp1 = MIN( jj+1, jpj )
1759	iim1 = MAX( ji-1, 1 )
1760	ijm1 = MAX( jj-1, 1 )
1761	IF( ( + bathy(iim1,ijp1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
1762	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
1763	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) THEN
1764	zenv(ji,jj) = rn_sbot_min
1765	ENDIF
1766	ENDIF
1767	END DO
1768	END DO
1769	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1770	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1771	!
1772	! smooth the bathymetry (if required)
1773	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1774	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1775	!
1776	jl = 0
1777	zrmax = 1._wp
1778	!
1779	!
1780	! set scaling factor used in reducing vertical gradients
1781	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1782	!
1783	! initialise temporary evelope depth arrays
1784	ztmpi1(:,:) = zenv(:,:)
1785	ztmpi2(:,:) = zenv(:,:)
1786	ztmpj1(:,:) = zenv(:,:)
1787	ztmpj2(:,:) = zenv(:,:)
1788	!
1789	! initialise temporary r-value arrays
1790	zri(:,:) = 1._wp
1791	zrj(:,:) = 1._wp
1792	! ! ================ !
1793	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1794	! ! ================ !
1795	jl = jl + 1
1796	zrmax = 0._wp
1797	! we set zrmax from previous r-values (zri and zrj) first
1798	! if set after current r-value calculation (as previously)
1799	! we could exit DO WHILE prematurely before checking r-value
1800	! of current zenv
1801	DO jj = 1, nlcj
1802	DO ji = 1, nlci
1803	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1804	END DO
1805	END DO
1806	zri(:,:) = 0._wp
1807	zrj(:,:) = 0._wp
1808	DO jj = 1, nlcj
1809	DO ji = 1, nlci
1810	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1811	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1812	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1813	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1814	END IF
1815	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1816	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1817	END IF
1818	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
1819	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
1820	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
1821	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
1822	END DO
1823	END DO
1824	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1825	!
1826	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
1827	!
1828	DO jj = 1, nlcj
1829	DO ji = 1, nlci
1830	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1831	END DO
1832	END DO
1833	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1834	CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1835	! ! ================ !
1836	END DO ! End loop !
1837	! ! ================ !
1838	DO jj = 1, jpj
1839	DO ji = 1, jpi
1840	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1841	END DO
1842	END DO
1843	!
1844	! Envelope bathymetry saved in hbatt
1845	hbatt(:,:) = zenv(:,:)
1846	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1847	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1848	DO jj = 1, jpj
1849	DO ji = 1, jpi
1850	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1851	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1852	END DO
1853	END DO
1854	ENDIF
1855	!
1856	IF(lwp) THEN ! Control print
1857	WRITE(numout,*)
1858	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1859	WRITE(numout,*)
1860	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1861	IF( nprint == 1 ) THEN
1862	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1863	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1864	ENDIF
1865	ENDIF
1866
1867	! ! ==============================
1868	! ! hbatu, hbatv, hbatf fields
1869	! ! ==============================
1870	IF(lwp) THEN
1871	WRITE(numout,*)
1872	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1873	ENDIF
1874	hbatu(:,:) = rn_sbot_min
1875	hbatv(:,:) = rn_sbot_min
1876	hbatf(:,:) = rn_sbot_min
1877	DO jj = 1, jpjm1
1878	DO ji = 1, jpim1 ! NO vector opt.
1879	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1880	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1881	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1882	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1883	END DO
1884	END DO
1885	!
1886	! Apply lateral boundary condition
1887	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1888	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
1889	DO jj = 1, jpj
1890	DO ji = 1, jpi
1891	IF( hbatu(ji,jj) == 0._wp ) THEN
1892	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1893	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1894	ENDIF
1895	END DO
1896	END DO
1897	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
1898	DO jj = 1, jpj
1899	DO ji = 1, jpi
1900	IF( hbatv(ji,jj) == 0._wp ) THEN
1901	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1902	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1903	ENDIF
1904	END DO
1905	END DO
1906	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
1907	DO jj = 1, jpj
1908	DO ji = 1, jpi
1909	IF( hbatf(ji,jj) == 0._wp ) THEN
1910	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1911	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1912	ENDIF
1913	END DO
1914	END DO
1915
1916	!!bug: key_helsinki a verifer
1917	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1918	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1919	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1920	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1921
1922	IF( nprint == 1 .AND. lwp ) THEN
1923	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1924	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1925	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1926	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1927	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1928	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1929	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1930	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1931	ENDIF
1932	!! helsinki
1933
1934	! ! =======================
1935	! ! s-ccordinate fields (gdep., e3.)
1936	! ! =======================
1937	!
1938	! non-dimensional "sigma" for model level depth at w- and t-levels
1939
1940
1941	!========================================================================
1942	! Song and Haidvogel 1994 (ln_s_sh94=T)
1943	! Siddorn and Furner 2012 (ln_sf12=T)
1944	! or tanh function (both false)
1945	!========================================================================
1946	IF ( ln_s_sh94 ) THEN
1947	CALL s_sh94()
1948	ELSE IF ( ln_s_sf12 ) THEN
1949	CALL s_sf12()
1950	ELSE
1951	CALL s_tanh()
1952	ENDIF
1953
1954	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
1955	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
1956	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
1957	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
1958	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
1959	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
1960	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1961
1962	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1963	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1964	!
1965	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
1966	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
1967	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
1968	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
1969	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
1970	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
1971	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
1972
1973	#if defined key_agrif
1974	! Ensure meaningful vertical scale factors in ghost lines/columns
1975	IF( .NOT. Agrif_Root() ) THEN
1976	!
1977	IF((nbondi == -1).OR.(nbondi == 2)) THEN
1978	e3u_0(1,:,:) = e3u_0(2,:,:)
1979	ENDIF
1980	!
1981	IF((nbondi == 1).OR.(nbondi == 2)) THEN
1982	e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
1983	ENDIF
1984	!
1985	IF((nbondj == -1).OR.(nbondj == 2)) THEN
1986	e3v_0(:,1,:) = e3v_0(:,2,:)
1987	ENDIF
1988	!
1989	IF((nbondj == 1).OR.(nbondj == 2)) THEN
1990	e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
1991	ENDIF
1992	!
1993	ENDIF
1994	#endif
1995
1996	fsdept(:,:,:) = gdept_0 (:,:,:)
1997	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1998	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1999	fse3t (:,:,:) = e3t_0 (:,:,:)
2000	fse3u (:,:,:) = e3u_0 (:,:,:)
2001	fse3v (:,:,:) = e3v_0 (:,:,:)
2002	fse3f (:,:,:) = e3f_0 (:,:,:)
2003	fse3w (:,:,:) = e3w_0 (:,:,:)
2004	fse3uw(:,:,:) = e3uw_0 (:,:,:)
2005	fse3vw(:,:,:) = e3vw_0 (:,:,:)
2006	!!
2007	! HYBRID :
2008	DO jj = 1, jpj
2009	DO ji = 1, jpi
2010	DO jk = 1, jpkm1
2011	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
2012	END DO
2013	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
2014	END DO
2015	END DO
2016	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
2017	& ' MAX ', MAXVAL( mbathy(:,:) )
2018
2019	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
2020	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
2021	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
2022	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
2023	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
2024	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
2025	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
2026	& ' w ', MINVAL( e3w_0 (:,:,:) )
2027
2028	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
2029	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
2030	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
2031	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
2032	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
2033	& ' w ', MAXVAL( e3w_0 (:,:,:) )
2034	ENDIF
2035	! END DO
2036	IF(lwp) THEN ! selected vertical profiles
2037	WRITE(numout,*)
2038	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2039	WRITE(numout,*) ' ~~~~~~ --------------------'
2040	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2041	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
2042	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2043	DO jj = mj0(20), mj1(20)
2044	DO ji = mi0(20), mi1(20)
2045	WRITE(numout,*)
2046	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2047	WRITE(numout,*) ' ~~~~~~ --------------------'
2048	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2049	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2050	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2051	END DO
2052	END DO
2053	DO jj = mj0(74), mj1(74)
2054	DO ji = mi0(100), mi1(100)
2055	WRITE(numout,*)
2056	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
2057	WRITE(numout,*) ' ~~~~~~ --------------------'
2058	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
2059	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
2060	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2061	END DO
2062	END DO
2063	ENDIF
2064
2065	!================================================================================
2066	! check the coordinate makes sense
2067	!================================================================================
2068	DO ji = 1, jpi
2069	DO jj = 1, jpj
2070
2071	IF( hbatt(ji,jj) > 0._wp) THEN
2072	DO jk = 1, mbathy(ji,jj)
2073	! check coordinate is monotonically increasing
2074	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
2075	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2076	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
2077	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
2078	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
2079	CALL ctl_stop( ctmp1 )
2080	ENDIF
2081	! and check it has never gone negative
2082	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
2083	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2084	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
2085	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2086	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2087	CALL ctl_stop( ctmp1 )
2088	ENDIF
2089	! and check it never exceeds the total depth
2090	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
2091	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2092	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
2093	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
2094	CALL ctl_stop( ctmp1 )
2095	ENDIF
2096	END DO
2097
2098	DO jk = 1, mbathy(ji,jj)-1
2099	! and check it never exceeds the total depth
2100	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
2101	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2102	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
2103	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
2104	CALL ctl_stop( ctmp1 )
2105	ENDIF
2106	END DO
2107
2108	ENDIF
2109
2110	END DO
2111	END DO
2112	!
2113	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2114	!
2115	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
2116	!
2117	END SUBROUTINE zgr_sco
2118
2119	!!======================================================================
2120	SUBROUTINE s_sh94()
2121
2122	!!----------------------------------------------------------------------
2123	!! * ROUTINE s_sh94 *
2124	!!
2125	!! ** Purpose : stretch the s-coordinate system
2126	!!
2127	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
2128	!! mixed S/sigma coordinate
2129	!!
2130	!! Reference : Song and Haidvogel 1994.
2131	!!----------------------------------------------------------------------
2132	!
2133	INTEGER :: ji, jj, jk ! dummy loop argument
2134	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2135	!
2136	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2137	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2138
2139	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2140	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2141
2142	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2143	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2144	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2145	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2146
2147	DO ji = 1, jpi
2148	DO jj = 1, jpj
2149
2150	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2151	DO jk = 1, jpk
2152	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2153	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2154	END DO
2155	ELSE ! shallow water, uniform sigma
2156	DO jk = 1, jpk
2157	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
2158	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2159	END DO
2160	ENDIF
2161	!
2162	DO jk = 1, jpkm1
2163	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2164	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2165	END DO
2166	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
2167	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2168	!
2169	! Coefficients for vertical depth as the sum of e3w scale factors
2170	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2171	DO jk = 2, jpk
2172	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2173	END DO
2174	!
2175	DO jk = 1, jpk
2176	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2177	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2178	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
2179	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
2180	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
2181	END DO
2182	!
2183	END DO ! for all jj's
2184	END DO ! for all ji's
2185
2186	DO ji = 1, jpim1
2187	DO jj = 1, jpjm1
2188	DO jk = 1, jpk
2189	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
2190	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2191	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
2192	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2193	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
2194	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
2195	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2196	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
2197	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
2198	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
2199	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
2200	!
2201	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2202	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2203	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2204	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2205	!
2206	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2207	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2208	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2209	END DO
2210	END DO
2211	END DO
2212
2213	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2214	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2215
2216	END SUBROUTINE s_sh94
2217
2218	SUBROUTINE s_sf12
2219
2220	!!----------------------------------------------------------------------
2221	!! * ROUTINE s_sf12 *
2222	!!
2223	!! ** Purpose : stretch the s-coordinate system
2224	!!
2225	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
2226	!! mixed S/sigma/Z coordinate
2227	!!
2228	!! This method allows the maintenance of fixed surface and or
2229	!! bottom cell resolutions (cf. geopotential coordinates)
2230	!! within an analytically derived stretched S-coordinate framework.
2231	!!
2232	!!
2233	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2234	!!----------------------------------------------------------------------
2235	!
2236	INTEGER :: ji, jj, jk ! dummy loop argument
2237	REAL(wp) :: zsmth ! smoothing around critical depth
2238	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
2239	!
2240	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2241	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
2242
2243	!
2244	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2245	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2246
2247	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
2248	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
2249	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
2250	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
2251
2252	DO ji = 1, jpi
2253	DO jj = 1, jpj
2254
2255	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2256
2257	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
2258	! could be changed by users but care must be taken to do so carefully
2259	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2260
2261	zzs = rn_zs / hbatt(ji,jj)
2262
2263	IF (rn_efold /= 0.0_wp) THEN
2264	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2265	ELSE
2266	zsmth = 1.0_wp
2267	ENDIF
2268
2269	DO jk = 1, jpk
2270	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2271	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2272	ENDDO
2273	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
2274	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
2275
2276	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2277
2278	DO jk = 1, jpk
2279	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
2280	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2281	END DO
2282
2283	ELSE ! shallow water, z coordinates
2284
2285	DO jk = 1, jpk
2286	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2287	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2288	END DO
2289
2290	ENDIF
2291
2292	DO jk = 1, jpkm1
2293	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2294	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2295	END DO
2296	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
2297	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2298
2299	! Coefficients for vertical depth as the sum of e3w scale factors
2300	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2301	DO jk = 2, jpk
2302	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2303	END DO
2304
2305	DO jk = 1, jpk
2306	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2307	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2308	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2309	END DO
2310
2311	ENDDO ! for all jj's
2312	ENDDO ! for all ji's
2313
2314	DO ji=1,jpi-1
2315	DO jj=1,jpj-1
2316
2317	DO jk = 1, jpk
2318	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
2319	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2320	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
2321	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2322	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
2323	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
2324	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2325	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
2326	( hbatt(ji,jj)+hbatt(ji+1,jj) )
2327	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
2328	( hbatt(ji,jj)+hbatt(ji,jj+1) )
2329
2330	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2331	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2332	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2333	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2334	!
2335	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2336	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2337	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2338	END DO
2339
2340	ENDDO
2341	ENDDO
2342	!
2343	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2344	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2345	CALL lbc_lnk(e3w_0 ,'T',1.)
2346	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2347	!
2348	! ! =============
2349
2350	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
2351	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2352
2353	END SUBROUTINE s_sf12
2354
2355	SUBROUTINE s_tanh()
2356
2357	!!----------------------------------------------------------------------
2358	!! * ROUTINE s_tanh*
2359	!!
2360	!! ** Purpose : stretch the s-coordinate system
2361	!!
2362	!! ** Method : s-coordinate stretch
2363	!!
2364	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2365	!!----------------------------------------------------------------------
2366
2367	INTEGER :: ji, jj, jk ! dummy loop argument
2368	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
2369
2370	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2371	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2372
2373	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2374	CALL wrk_alloc( jpk, z_esigt, z_esigw )
2375
2376	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
2377	z_esigt = 0._wp ; z_esigw = 0._wp
2378
2379	DO jk = 1, jpk
2380	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2381	z_gsigt(jk) = -fssig( REAL(jk,wp) )
2382	END DO
2383	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
2384	!
2385	! Coefficients for vertical scale factors at w-, t- levels
2386	!!gm bug : define it from analytical function, not like juste bellow....
2387	!!gm or betteroffer the 2 possibilities....
2388	DO jk = 1, jpkm1
2389	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
2390	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2391	END DO
2392	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
2393	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2394	!
2395	! Coefficients for vertical depth as the sum of e3w scale factors
2396	z_gsi3w(1) = 0.5_wp * z_esigw(1)
2397	DO jk = 2, jpk
2398	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2399	END DO
2400	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2401	DO jk = 1, jpk
2402	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2403	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2404	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
2405	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
2406	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
2407	END DO
2408	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2409	DO jj = 1, jpj
2410	DO ji = 1, jpi
2411	DO jk = 1, jpk
2412	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2413	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2414	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2415	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2416	!
2417	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2418	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2419	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2420	END DO
2421	END DO
2422	END DO
2423
2424	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
2425	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
2426
2427	END SUBROUTINE s_tanh
2428
2429	FUNCTION fssig( pk ) RESULT( pf )
2430	!!----------------------------------------------------------------------
2431	!! * ROUTINE fssig *
2432	!!
2433	!! ** Purpose : provide the analytical function in s-coordinate
2434	!!
2435	!! ** Method : the function provide 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) :: pk ! continuous "k" coordinate
2441	REAL(wp) :: pf ! sigma value
2442	!!----------------------------------------------------------------------
2443	!
2444	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2445	& - TANH( rn_thetb * rn_theta ) ) &
2446	& * ( COSH( rn_theta ) &
2447	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2448	& / ( 2._wp * SINH( rn_theta ) )
2449	!
2450	END FUNCTION fssig
2451
2452
2453	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2454	!!----------------------------------------------------------------------
2455	!! * ROUTINE fssig1 *
2456	!!
2457	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2458	!!
2459	!! ** Method : the function provides the non-dimensional position of
2460	!! T and W (i.e. between 0 and 1)
2461	!! T-points at integer values (between 1 and jpk)
2462	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2463	!!----------------------------------------------------------------------
2464	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2465	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2466	REAL(wp) :: pf1 ! sigma value
2467	!!----------------------------------------------------------------------
2468	!
2469	IF ( rn_theta == 0 ) then ! uniform sigma
2470	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2471	ELSE ! stretched sigma
2472	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2473	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2474	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2475	ENDIF
2476	!
2477	END FUNCTION fssig1
2478
2479
2480	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2481	!!----------------------------------------------------------------------
2482	!! * ROUTINE fgamma *
2483	!!
2484	!! ** Purpose : provide analytical function for the s-coordinate
2485	!!
2486	!! ** Method : the function provides the non-dimensional position of
2487	!! T and W (i.e. between 0 and 1)
2488	!! T-points at integer values (between 1 and jpk)
2489	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2490	!!
2491	!! This method allows the maintenance of fixed surface and or
2492	!! bottom cell resolutions (cf. geopotential coordinates)
2493	!! within an analytically derived stretched S-coordinate framework.
2494	!!
2495	!! Reference : Siddorn and Furner, in prep
2496	!!----------------------------------------------------------------------
2497	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2498	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2499	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2500	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2501	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2502	REAL(wp) :: za1,za2,za3 ! local variables
2503	REAL(wp) :: zn1,zn2 ! local variables
2504	REAL(wp) :: za,zb,zx ! local variables
2505	integer :: jk
2506	!!----------------------------------------------------------------------
2507	!
2508
2509	zn1 = 1./(jpk-1.)
2510	zn2 = 1. - zn1
2511
2512	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2513	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2514	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2515
2516	za = pzb - za3*(pzs-za1)-za2
2517	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2518	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2519	zx = 1.0_wp-za/2.0_wp-zb
2520
2521	DO jk = 1, jpk
2522	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2523	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2524	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2525	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2526	ENDDO
2527
2528	!
2529	END FUNCTION fgamma
2530
2531	!!======================================================================
2532	END MODULE domzgr

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

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

Context Navigation

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

Download in other formats: