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

Context Navigation

source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5315

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

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

Download in other formats: