Context Navigation

domzgr.F90 @ 7452

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

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

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

Context Navigation

source: branches/UKMO/r6232_CO6_CO5_zenv_pomsdwl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7452

Download in other formats: