Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domzgr.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 136.4 KB

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

Download in other formats: