Context Navigation

domzgr.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

Property svn:keywords set to Id

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

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

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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4990

Download in other formats: