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

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

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

Context Navigation

source: branches/UKMO/2015_V36_STABLE_CO6_CO5_zenv_pomsdwl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5654

Last change on this file since 5654 was 5654, checked in by deazer, 9 years ago
Added in changes to allow CO5 like run to verify CO6 with NEMO STABLE VN3.6
File size: 130.3 KB

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

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

Download in other formats: