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/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC – NEMO

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/domzgr.F90 @ 6895

Last change on this file since 6895 was 6895, checked in by flavoni, 8 years ago
commit usrdef routines for OVERFLOW
File size: 113.3 KB

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

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

Download in other formats: