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

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

domzgr.F90 in branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6012

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

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

Download in other formats: