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

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

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

Context Navigation

source: branches/UKMO/dev_r5107_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5473

Last change on this file since 5473 was 5473, checked in by cguiavarch, 9 years ago
Clear svn keywords from UKMO/dev_r5107_hadgem3_cplfld
File size: 127.9 KB

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

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

Download in other formats: