Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domzgr.F90 @ 6510

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6510

Download in other formats: