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.

grid_zgr.f90 in NEMO/trunk/tools/SIREN/src – NEMO

Context Navigation

source: NEMO/trunk/tools/SIREN/src/grid_zgr.f90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago
Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing Standardisation of routines header with version 4.0 and year 2018 Fix for some broken symlinks
File size: 192.3 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! MODULE: grid_zgr
6	!
7	! DESCRIPTION:
8	!> @brief This module manage Vertical grid.
9	!>
10	!> @details
11	!> ** Purpose : set the depth of model levels and the resulting
12	!> vertical scale factors.
13	!>
14	!> ** Method :
15	!> - reference 1D vertical coordinate (gdep._1d, e3._1d)
16	!> - read/set ocean depth and ocean levels (bathy, mbathy)
17	!> - vertical coordinate (gdep., e3.) depending on the coordinate chosen :
18	!> - ln_zco=T z-coordinate
19	!> - ln_zps=T z-coordinate with partial steps
20	!> - ln_zco=T s-coordinate
21	!>
22	!> ** Action : define gdep., e3., mbathy and bathy
23	!>
24	!> @author
25	!> G, Madec
26	! REVISION HISTORY:
27	!> @date December, 1995 - Original code : s vertical coordinate
28	!> @date July, 1997
29	!> - lbc_lnk call
30	!> @date September, 2002
31	!> - A. Bozec, G. Madec : F90: Free form and module
32	!> @date September, 2002
33	!> - A. de Miranda : rigid-lid + islands
34	!> @date August, 2003
35	!> - G. Madec : Free form and module
36	!> @date October, 2005
37	!> - A. Beckmann : modifications for hybrid s-ccordinates & new stretching function
38	!> @date April, 2006
39	!> - R. Benshila, G. Madec : add zgr_zco
40	!> @date June, 2008
41	!> - G. Madec : insertion of domzgr_zps.h90 & conding style
42	!> @date July, 2009
43	!> - R. Benshila : Suppression of rigid-lid option
44	!> @date November, 2011
45	!> - G. Madec : add mbk. arrays associated to the deepest ocean level
46	!> @date August, 2012
47	!> - J. Siddorn : added Siddorn and Furner stretching function
48	!> @date December, 2012
49	!> - R. Bourdalle-Badie and G. Reffray : modify C1D case
50	!> @date November, 2014
51	!> - P. Mathiot and C. Harris : add ice shelf capabilitye
52	!> @date November, 2015
53	!> - H. Liu : Modifications for Wetting/Drying
54	!> @date October, 2016
55	!> - J, Paul : update from trunk (revision 6961): add wetting and drying, ice sheet coupling..
56	!> - J, Paul : do not use anymore special case for ORCA grid.
57	!> @date November, 2016
58	!> - J, Paul : vertical scale factors e3. = dk[gdep] or old definition
59	!>
60	!> @note Software governed by the CeCILL licence (./LICENSE)
61	!----------------------------------------------------------------------
62	MODULE grid_zgr
63	USE netcdf ! nf90 library
64	USE kind ! F90 kind parameter
65	USE fct ! basic usefull function
66	USE global ! global parameter
67	USE phycst ! physical constant
68	USE logger ! log file manager
69	USE file ! file manager
70	USE var ! variable manager
71	USE dim ! dimension manager
72	USE dom ! domain manager
73	USE grid ! grid manager
74	USE iom ! I/O manager
75	USE mpp ! MPP manager
76	USE iom_mpp ! I/O MPP manager
77	USE lbc ! lateral boundary conditions
78	USE grid_hgr ! Horizontal mesh
79	IMPLICIT NONE
80	! NOTE_avoid_public_variables_if_possible
81	! type and variable
82	PUBLIC :: TNAMZ
83
84	PUBLIC :: tg_gdepw_1d
85	PUBLIC :: tg_gdept_1d
86	PUBLIC :: tg_e3w_1d
87	PUBLIC :: tg_e3t_1d
88	PUBLIC :: tg_e3tp
89	PUBLIC :: tg_e3wp
90
91	PUBLIC :: tg_rx1
92
93	PUBLIC :: tg_mbathy
94	PUBLIC :: tg_misfdep
95
96	PUBLIC :: tg_gdept_0
97	PUBLIC :: tg_gdepw_0
98	! PUBLIC :: tg_gdep3w_0 !useless to create meshmask
99	PUBLIC :: tg_e3t_0
100	PUBLIC :: tg_e3u_0
101	PUBLIC :: tg_e3v_0
102	PUBLIC :: tg_e3w_0
103	PUBLIC :: tg_e3f_0 !useless to create meshmask
104	PUBLIC :: tg_e3uw_0 !useless to create meshmask
105	PUBLIC :: tg_e3vw_0 !useless to create meshmask
106
107	PUBLIC :: tg_mbkt
108	! PUBLIC :: tg_mbku !useless to create meshmask
109	! PUBLIC :: tg_mbkv !useless to create meshmask
110	PUBLIC :: tg_mikt
111	! PUBLIC :: tg_miku !useless to create meshmask
112	! PUBLIC :: tg_mikv !useless to create meshmask
113	! PUBLIC :: tg_mikf !useless to create meshmask
114
115	PUBLIC :: tg_hbatt ! sco
116	PUBLIC :: tg_hbatu ! sco
117	PUBLIC :: tg_hbatv ! sco
118	PUBLIC :: tg_hbatf ! sco
119
120	PUBLIC :: tg_gsigt ! sco(tanh)
121	PUBLIC :: tg_gsigw ! sco(tanh)
122	PUBLIC :: tg_gsi3w ! sco(tanh)
123	PUBLIC :: tg_esigt ! sco(tanh)
124	PUBLIC :: tg_esigw ! sco(tanh)
125
126	! function and subroutine
127	PUBLIC :: grid_zgr_init
128	PUBLIC :: grid_zgr_nam
129	PUBLIC :: grid_zgr_fill
130	PUBLIC :: grid_zgr_clean
131
132	PUBLIC :: grid_zgr_zps_init
133	PUBLIC :: grid_zgr_zps_clean
134	PUBLIC :: grid_zgr_sco_init
135	PUBLIC :: grid_zgr_sco_clean
136	PUBLIC :: grid_zgr_sco_stiff
137
138	PRIVATE :: grid_zgr__z
139	PRIVATE :: grid_zgr__bat
140	PRIVATE :: grid_zgr__zco
141	! PRIVATE :: grid_zgr__bat_zoom
142	PRIVATE :: grid_zgr__bat_ctl
143	PRIVATE :: grid_zgr__bot_level
144	PRIVATE :: grid_zgr__top_level
145	PRIVATE :: grid_zgr__zps_fill
146	PRIVATE :: grid_zgr__isf_fill
147	! PRIVATE :: grid_zgr__isf_fill_e3x
148	PRIVATE :: grid_zgr__isf_fill_e3uw
149	! PRIVATE :: grid_zgr__isf_fill_gdep3w_0
150	PRIVATE :: grid_zgr__sco_fill
151	PRIVATE :: grid_zgr__sco_s_sh94
152	PRIVATE :: grid_zgr__sco_s_sf12
153	PRIVATE :: grid_zgr__sco_s_tanh
154	PRIVATE :: grid_zgr__sco_fssig !: tanh stretch function
155	PRIVATE :: grid_zgr__sco_fssig1 !: Song and Haidvogel 1994 stretch function
156	PRIVATE :: grid_zgr__sco_fgamma !: Siddorn and Furner 2012 stretching function
157
158	TYPE TNAMZ
159
160	CHARACTER(LEN=lc) :: c_coord
161	INTEGER(i4) :: i_perio
162
163	LOGICAL :: l_zco
164	LOGICAL :: l_zps
165	LOGICAL :: l_sco
166	LOGICAL :: l_isfcav
167	LOGICAL :: l_iscpl
168	LOGICAL :: l_wd
169	INTEGER(i4) :: i_nlevel
170
171	REAL(dp) :: d_ppsur
172	REAL(dp) :: d_ppa0
173	REAL(dp) :: d_ppa1
174	REAL(dp) :: d_ppkth
175	REAL(dp) :: d_ppacr
176	REAL(dp) :: d_ppdzmin
177	REAL(dp) :: d_pphmax
178	LOGICAL :: l_dbletanh
179	REAL(dp) :: d_ppa2
180	REAL(dp) :: d_ppkth2
181	REAL(dp) :: d_ppacr2
182
183	REAL(dp) :: d_hmin
184	REAL(dp) :: d_isfhmin
185
186	REAL(dp) :: d_e3zps_min
187	REAL(dp) :: d_e3zps_rat
188	! INTEGER(i4) :: i_msh
189
190	LOGICAL :: l_s_sh94
191	LOGICAL :: l_s_sf12
192	REAL(dp) :: d_sbot_min
193	REAL(dp) :: d_sbot_max
194	! Song and Haidvogel 1994 stretching additional parameters
195	REAL(dp) :: d_rmax
196	REAL(dp) :: d_hc
197	REAL(dp) :: d_theta
198	REAL(dp) :: d_thetb
199	REAL(dp) :: d_bb
200	! Siddorn and Furner stretching additional parameters
201	LOGICAL :: l_sigcrit
202	REAL(dp) :: d_alpha
203	REAL(dp) :: d_efold
204	REAL(dp) :: d_zs
205	REAL(dp) :: d_zb_a
206	REAL(dp) :: d_zb_b
207
208	INTEGER(i4) :: i_cla
209
210	REAL(dp) :: d_wdmin1
211	REAL(dp) :: d_wdmin2
212	REAL(dp) :: d_wdld
213
214	! CHARACTER(LEN=lc) :: c_cfg
215	! INTEGER(i4) :: i_cfg
216	! INTEGER(i4) :: i_bench
217	! LOGICAL :: l_zoom
218	LOGICAL :: l_c1d
219	LOGICAL :: l_e3_dep
220
221	! CHARACTER(LEN=lc) :: c_cfz
222	! INTEGER(i4) :: i_izoom
223	! INTEGER(i4) :: i_jzoom
224	! LOGICAL :: l_zoom_s
225	! LOGICAL :: l_zoom_e
226	! LOGICAL :: l_zoom_w
227	! LOGICAL :: l_zoom_n
228
229	END TYPE
230
231	TYPE(TVAR), SAVE :: tg_gdepw_1d !zco & zps & sco
232	TYPE(TVAR), SAVE :: tg_gdept_1d !zco & zps & sco
233	TYPE(TVAR), SAVE :: tg_e3w_1d !zco & zps
234	TYPE(TVAR), SAVE :: tg_e3t_1d !zco & zps
235	TYPE(TVAR), SAVE :: tg_e3tp ! zps
236	TYPE(TVAR), SAVE :: tg_e3wp ! zps
237
238	TYPE(TVAR), SAVE :: tg_rx1 ! sco
239
240	TYPE(TVAR), SAVE :: tg_mbathy !zco & zps & sco
241	TYPE(TVAR), SAVE :: tg_misfdep
242
243	TYPE(TVAR), SAVE :: tg_gdept_0 ! zps & sco
244	TYPE(TVAR), SAVE :: tg_gdepw_0 ! zps & sco
245	!TYPE(TVAR), SAVE :: tg_gdep3w_0
246	TYPE(TVAR), SAVE :: tg_e3t_0 ! zps & sco
247	TYPE(TVAR), SAVE :: tg_e3u_0 ! zps & sco
248	TYPE(TVAR), SAVE :: tg_e3v_0 ! zps & sco
249	TYPE(TVAR), SAVE :: tg_e3w_0 ! zps & sco
250	TYPE(TVAR), SAVE :: tg_e3f_0
251	TYPE(TVAR), SAVE :: tg_e3uw_0
252	TYPE(TVAR), SAVE :: tg_e3vw_0
253
254	TYPE(TVAR), SAVE :: tg_mbkt !zco & zps & sco
255	!TYPE(TVAR), SAVE :: tg_mbku
256	!TYPE(TVAR), SAVE :: tg_mbkv
257	TYPE(TVAR), SAVE :: tg_mikt !zco & zps & sco
258	!TYPE(TVAR), SAVE :: tg_miku
259	!TYPE(TVAR), SAVE :: tg_mikv
260	!TYPE(TVAR), SAVE :: tg_mikf
261
262	TYPE(TVAR), SAVE :: tg_hbatt ! sco
263	TYPE(TVAR), SAVE :: tg_hbatu ! sco
264	TYPE(TVAR), SAVE :: tg_hbatv ! sco
265	TYPE(TVAR), SAVE :: tg_hbatf ! sco
266
267	TYPE(TVAR), SAVE :: tg_gsigt ! sco(tanh)
268	TYPE(TVAR), SAVE :: tg_gsigw ! sco(tanh)
269	TYPE(TVAR), SAVE :: tg_gsi3w ! sco(tanh)
270	TYPE(TVAR), SAVE :: tg_esigt ! sco(tanh)
271	TYPE(TVAR), SAVE :: tg_esigw ! sco(tanh)
272
273	CONTAINS
274	!-------------------------------------------------------------------
275	!> @brief This function initialise global variable needed to compute vertical
276	!> mesh
277	!>
278	!> @author J.Paul
279	!> @date September, 2015 - Initial version
280	!>
281	!> @param[in] jpi
282	!> @param[in] jpj
283	!> @param[in] jpk
284	!> @param[in] ld_sco
285	!-------------------------------------------------------------------
286	SUBROUTINE grid_zgr_init( jpi,jpj,jpk, ld_sco )
287	IMPLICIT NONE
288	! Argument
289	INTEGER(i4), INTENT(IN) :: jpi
290	INTEGER(i4), INTENT(IN) :: jpj
291	INTEGER(i4), INTENT(IN) :: jpk
292	LOGICAL , INTENT(IN) :: ld_sco
293
294	! local variable
295	REAL(dp), DIMENSION(jpk) :: dl_tmp1D
296	REAL(dp), DIMENSION(jpi,jpj) :: dl_tmp2D
297	REAL(dp), DIMENSION(jpi,jpj,jpk) :: dl_tmp3D
298	! loop indices
299	!----------------------------------------------------------------
300
301	! variable 1D
302	dl_tmp1D(:) =dp_fill
303
304	tg_gdepw_1d=var_init('gdepw_1d',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
305	tg_gdept_1d=var_init('gdept_1d',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
306	tg_e3w_1d =var_init('e3w_1d ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
307	tg_e3t_1d =var_init('e3t_1d ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
308
309	!only sco
310	IF( ld_sco )THEN
311	tg_gsigt =var_init('gsigt ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
312	tg_gsigw =var_init('gsigw ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
313	tg_gsi3w =var_init('gsi3w ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
314	tg_esigt =var_init('esigt ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
315	tg_esigw =var_init('esigw ',dl_tmp1D(:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
316	ENDIF
317
318	! variable 2D
319	dl_tmp2D(:,:) =dp_fill_i2
320
321	tg_mbkt =var_init('mbkt ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
322	!tg_mbku =var_init('mbku ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
323	!tg_mbkv =var_init('mbkv ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
324	tg_mikt =var_init('mikt ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
325	!tg_miku =var_init('miku ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
326	!tg_mikv =var_init('mikv ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
327	!tg_mikf =var_init('mikf ',dl_tmp2D(:,:) , dd_fill=dp_fill_i2, id_type=NF90_SHORT)
328
329	dl_tmp2D(:,:) =dp_fill_i4
330
331	tg_mbathy =var_init('mbathy ',dl_tmp2D(:,:) , dd_fill=dp_fill_i4, id_type=NF90_INT)
332	tg_misfdep =var_init('misfdep ',dl_tmp2D(:,:) , dd_fill=dp_fill_i4, id_type=NF90_INT)
333
334	dl_tmp2D(:,:) =dp_fill
335
336	! only sco
337	IF( ld_sco )THEN
338	tg_hbatt =var_init('hbatt ',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
339	tg_hbatu =var_init('hbatu ',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
340	tg_hbatv =var_init('hbatv ',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
341	tg_hbatf =var_init('hbatf ',dl_tmp2D(:,:) , dd_fill=dp_fill, id_type=NF90_DOUBLE)
342	ENDIF
343
344	! variable 3D
345	dl_tmp3D(:,:,:)=dp_fill
346
347	tg_gdept_0 =var_init('gdept_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
348	tg_gdepw_0 =var_init('gdepw_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
349	!tg_gdep3w_0=var_init('gdep3w_0',dl_tmp3D(:,:,:), dd_fill=dp_fill_sp, id_type=NF90_FLOAT)
350
351	tg_e3t_0 =var_init('e3t_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
352	tg_e3u_0 =var_init('e3u_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
353	tg_e3v_0 =var_init('e3v_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
354	tg_e3w_0 =var_init('e3w_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
355	tg_e3f_0 =var_init('e3f_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
356	tg_e3uw_0 =var_init('e3uw_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
357	tg_e3vw_0 =var_init('e3vw_0 ',dl_tmp3D(:,:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
358
359	END SUBROUTINE grid_zgr_init
360	!-------------------------------------------------------------------
361	!> @brief This function clean hgr structure
362	!>
363	!> @author J.Paul
364	!> @date September, 2015 - Initial version
365	!>
366	!> @param[in] ld_sco
367	!-------------------------------------------------------------------
368	SUBROUTINE grid_zgr_clean(ld_sco)
369	IMPLICIT NONE
370	! Argument
371	LOGICAL , INTENT(IN) :: ld_sco
372
373	! local variable
374	! loop indices
375	!----------------------------------------------------------------
376
377	CALL var_clean(tg_gdepw_1d)
378	CALL var_clean(tg_gdept_1d)
379	CALL var_clean(tg_e3w_1d )
380	CALL var_clean(tg_e3t_1d )
381
382	IF( ld_sco )THEN
383	CALL var_clean(tg_gsigt )
384	CALL var_clean(tg_gsigw )
385	CALL var_clean(tg_gsi3w )
386	CALL var_clean(tg_esigt )
387	CALL var_clean(tg_esigw )
388	ENDIF
389
390	CALL var_clean(tg_mbathy )
391	CALL var_clean(tg_misfdep )
392
393	CALL var_clean(tg_mbkt )
394	!CALL var_clean(tg_mbku )
395	!CALL var_clean(tg_mbkv )
396	CALL var_clean(tg_mikt )
397	!CALL var_clean(tg_miku )
398	!CALL var_clean(tg_mikv )
399	!CALL var_clean(tg_mikf )
400
401	IF( ld_sco )THEN
402	CALL var_clean(tg_hbatt )
403	CALL var_clean(tg_hbatu )
404	CALL var_clean(tg_hbatv )
405	CALL var_clean(tg_hbatf )
406	ENDIF
407
408	CALL var_clean(tg_gdept_0 )
409	CALL var_clean(tg_gdepw_0 )
410	!CALL var_clean(tg_gdep3w_0)
411
412	CALL var_clean(tg_e3t_0 )
413	CALL var_clean(tg_e3u_0 )
414	CALL var_clean(tg_e3v_0 )
415	CALL var_clean(tg_e3f_0 )
416	CALL var_clean(tg_e3uw_0 )
417	CALL var_clean(tg_e3vw_0 )
418
419	END SUBROUTINE grid_zgr_clean
420	!-------------------------------------------------------------------
421	!> @brief This function initialise zgr namelist structure
422	!>
423	!> @author J.Paul
424	!> @date September, 2015 - Initial version
425	!>
426	!> @param[in] cd_coord
427	!> @param[in] id_perio
428	!> @param[in] cd_namelist
429	!> @return hgr namelist structure
430	!-------------------------------------------------------------------
431	FUNCTION grid_zgr_nam( cd_coord,id_perio,cd_namelist )
432	IMPLICIT NONE
433	! Argument
434	CHARACTER(LEN=*), INTENT(IN) :: cd_coord
435	INTEGER(i4) , INTENT(IN) :: id_perio
436
437	CHARACTER(LEN=*), INTENT(IN) :: cd_namelist
438
439	! function
440	TYPE(TNAMZ) :: grid_zgr_nam
441
442	! local variable
443	INTEGER(i4) :: il_status
444	INTEGER(i4) :: il_fileid
445
446	LOGICAL :: ll_exist
447
448	! namelist
449	! namzgr
450	LOGICAL :: ln_zco = .FALSE.
451	LOGICAL :: ln_zps = .FALSE.
452	LOGICAL :: ln_sco = .FALSE.
453	LOGICAL :: ln_isfcav = .FALSE.
454	LOGICAL :: ln_iscpl = .FALSE.
455	LOGICAL :: ln_wd = .FALSE.
456	INTEGER(i4) :: in_nlevel = 75
457
458	! namdmin
459	REAL(dp) :: dn_hmin = NF90_FILL_DOUBLE
460	REAL(dp) :: dn_isfhmin = NF90_FILL_DOUBLE
461
462	! namzco
463	REAL(dp) :: dn_ppsur = NF90_FILL_DOUBLE
464	REAL(dp) :: dn_ppa0 = NF90_FILL_DOUBLE
465	REAL(dp) :: dn_ppa1 = NF90_FILL_DOUBLE
466	REAL(dp) :: dn_ppkth = NF90_FILL_DOUBLE
467	REAL(dp) :: dn_ppacr = NF90_FILL_DOUBLE
468	REAL(dp) :: dn_ppdzmin = NF90_FILL_DOUBLE
469	REAL(dp) :: dn_pphmax = NF90_FILL_DOUBLE
470	LOGICAL :: ln_dbletanh = .FALSE.
471	REAL(dp) :: dn_ppa2 = NF90_FILL_DOUBLE
472	REAL(dp) :: dn_ppkth2 = NF90_FILL_DOUBLE
473	REAL(dp) :: dn_ppacr2 = NF90_FILL_DOUBLE
474
475	! namzps
476	REAL(dp) :: dn_e3zps_min= NF90_FILL_DOUBLE
477	REAL(dp) :: dn_e3zps_rat= NF90_FILL_DOUBLE
478	! INTEGER(i4) :: in_msh = NF90_FILL_INT
479
480	! namsco
481	LOGICAL :: ln_s_sh94 = .FALSE.
482	LOGICAL :: ln_s_sf12 = .FALSE.
483	REAL(dp) :: dn_sbot_min = NF90_FILL_DOUBLE
484	REAL(dp) :: dn_sbot_max = NF90_FILL_DOUBLE
485	REAL(dp) :: dn_rmax = NF90_FILL_DOUBLE
486	REAL(dp) :: dn_hc = NF90_FILL_DOUBLE
487	!
488	REAL(dp) :: dn_theta = NF90_FILL_DOUBLE
489	REAL(dp) :: dn_thetb = NF90_FILL_DOUBLE
490	REAL(dp) :: dn_bb = NF90_FILL_DOUBLE
491	!
492	LOGICAL :: ln_sigcrit = .FALSE.
493	REAL(dp) :: dn_alpha = NF90_FILL_DOUBLE
494	REAL(dp) :: dn_efold = NF90_FILL_DOUBLE
495	REAL(dp) :: dn_zs = NF90_FILL_DOUBLE
496	REAL(dp) :: dn_zb_a = NF90_FILL_DOUBLE
497	REAL(dp) :: dn_zb_b = NF90_FILL_DOUBLE
498
499	! ! namcla
500	! INTEGER(i4) :: in_cla = 0
501
502	! namwd
503	REAL(dp) :: dn_wdmin1 = NF90_FILL_DOUBLE
504	REAL(dp) :: dn_wdmin2 = NF90_FILL_DOUBLE
505	REAL(dp) :: dn_wdld = NF90_FILL_DOUBLE
506
507	! namgrd
508	! CHARACTER(LEN=lc) :: cn_cfg = ''
509	! INTEGER(i4) :: in_cfg = 0
510	! INTEGER(i4) :: in_bench = 0
511	! LOGICAL :: ln_zoom = .FALSE.
512	LOGICAL :: ln_c1d = .FALSE.
513	LOGICAL :: ln_e3_dep = .FALSE.
514
515	! ! namzoom
516	! CHARACTER(LEN=lc) :: cn_cfz =''
517	! INTEGER(i4) :: in_izoom = NF90_FILL_INT
518	! INTEGER(i4) :: in_jzoom = NF90_FILL_INT
519	! LOGICAL :: ln_zoom_s = .FALSE.
520	! LOGICAL :: ln_zoom_e = .FALSE.
521	! LOGICAL :: ln_zoom_w = .FALSE.
522	! LOGICAL :: ln_zoom_n = .FALSE.
523	!----------------------------------------------------------------
524	NAMELIST /namzgr/ &
525	& ln_zco, & !< z-coordinate
526	& ln_zps, & !< z-coordinate with partial steps
527	& ln_sco, & !< s-coordinate
528	& ln_isfcav, & !< presence of ISF
529	& ln_iscpl, & !< coupling with ice sheet
530	& ln_wd, & !< Wetting/drying activation
531	& in_nlevel !< number of vertical level
532
533	NAMELIST /namdmin/ &
534	& dn_hmin, & !< minimum ocean depth (>0) or minimum number of ocean levels (<0)
535	& dn_isfhmin !< threshold to discriminate grounded ice to floating ice
536
537	NAMELIST /namzco/ &
538	& dn_ppsur, &
539	& dn_ppa0, &
540	& dn_ppa1, &
541	& dn_ppkth, &
542	& dn_ppacr, &
543	& dn_ppdzmin, &
544	& dn_pphmax, &
545	& ln_dbletanh, &
546	& dn_ppa2, &
547	& dn_ppkth2, &
548	& dn_ppacr2
549
550	NAMELIST /namzps/ &
551	& dn_e3zps_min, &
552	& dn_e3zps_rat!, &
553	! & in_msh
554
555	NAMELIST /namsco/ &
556	& ln_s_sh94, & !< use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1
557	& ln_s_sf12, & !< use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma
558	& dn_sbot_min, & !< minimum depth of s-bottom surface (>0) (m)
559	& dn_sbot_max, & !< maximum depth of s-bottom surface (= ocean depth) (>0) (m)
560	& dn_hc, & !< Critical depth for transition from sigma to stretched coordinates
561	! Song and Haidvogel 1994 stretching parameters
562	& dn_rmax, & !< maximum cut-off r-value allowed (0<dn_rmax<1)
563	& dn_theta, & !< surface control parameter (0<=dn_theta<=20)
564	& dn_thetb, & !< bottom control parameter (0<=dn_thetb<= 1)
565	& dn_bb, & !< stretching parameter ( dn_bb=0; top only, dn_bb =1; top and bottom)
566	! Siddorn and Furner stretching parameters
567	& ln_sigcrit, & !< use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
568	& dn_alpha, & !< control parameter ( > 1 stretch towards surface, < 1 towards seabed)
569	& dn_efold, & !< efold length scale for transition to stretched coord
570	& dn_zs, & !< depth of surface grid box
571	!< bottom cell depth (Zb) is a linear function of water depth Zb = H*rn_zb_a + rn_zb_b'
572	& dn_zb_a, & !< bathymetry scaling factor for calculating Zb
573	& dn_zb_b !< offset for calculating Zb
574
575	! NAMELIST /namcla/ &
576	! & in_cla !< =1 cross land advection for exchanges through some straits (ORCA2)
577
578	NAMELIST /namwd/ & !< wetting and drying
579	& dn_wdmin1, & !< minimum water depth on dried cells
580	& dn_wdmin2, & !< tolerrance of minimum water depth on dried cells
581	& dn_wdld !< land elevation below which wetting/drying
582
583	NAMELIST/namgrd/ & !< orca grid namelist
584	! & cn_cfg, & !< name of the configuration (orca)
585	! & in_cfg, & !< resolution of the configuration (2,1,025..)
586	! & in_bench, & !< benchmark parameter (in_mshhgr = 5 )
587	! & ln_zoom, & !< use zoom
588	& ln_c1d, & !< use configuration 1D
589	& ln_e3_dep !< new vertical scale factors [T, F:old definition]
590
591	! NAMELIST /namzoom/&
592	! & cn_cfz, & !< name of the zoom of configuration
593	! & in_izoom, & !< left bottom i-indices of the zoom in data domain indices
594	! & in_jzoom, & !< left bottom j-indices of the zoom in data domain indices
595	! & ln_zoom_s, & !< South zoom type flag
596	! & ln_zoom_e, & !< East zoom type flag
597	! & ln_zoom_w, & !< West zoom type flag
598	! & ln_zoom_n !< North zoom type flag
599	!----------------------------------------------------------------
600	!1-2 read namelist
601	INQUIRE(FILE=TRIM(cd_namelist), EXIST=ll_exist)
602	IF( ll_exist )THEN
603
604	il_fileid=fct_getunit()
605
606	OPEN( il_fileid, FILE=TRIM(cd_namelist), &
607	& FORM='FORMATTED', &
608	& ACCESS='SEQUENTIAL', &
609	& STATUS='OLD', &
610	& ACTION='READ', &
611	& IOSTAT=il_status)
612	CALL fct_err(il_status)
613	IF( il_status /= 0 )THEN
614	CALL logger_fatal("GRID ZGR NAM: error opening "//&
615	& TRIM(cd_namelist))
616	ENDIF
617
618	READ( il_fileid, NML = namzgr )
619	READ( il_fileid, NML = namdmin )
620	READ( il_fileid, NML = namzco )
621
622	IF( ln_zps ) READ( il_fileid, NML = namzps )
623	IF( ln_sco ) READ( il_fileid, NML = namsco )
624	! READ( il_fileid, NML = namcla )
625	READ( il_fileid, NML = namwd )
626	READ( il_fileid, NML = namgrd )
627	! IF( ln_zoom ) READ( il_fileid, NML = namzoom )
628
629	CLOSE( il_fileid, IOSTAT=il_status )
630	CALL fct_err(il_status)
631	IF( il_status /= 0 )THEN
632	CALL logger_error("GRID ZGR NAM: closing "//TRIM(cd_namelist))
633	ENDIF
634
635	grid_zgr_nam%c_coord = TRIM(cd_coord)
636	grid_zgr_nam%i_perio = id_perio
637
638	grid_zgr_nam%l_zco = ln_zco
639	grid_zgr_nam%l_zps = ln_zps
640	grid_zgr_nam%l_sco = ln_sco
641	grid_zgr_nam%l_isfcav = ln_isfcav
642	grid_zgr_nam%l_iscpl = ln_iscpl
643	grid_zgr_nam%l_wd = ln_wd
644	grid_zgr_nam%i_nlevel = in_nlevel
645
646	grid_zgr_nam%d_hmin = dn_hmin
647	grid_zgr_nam%d_isfhmin = dn_isfhmin
648
649	grid_zgr_nam%d_ppsur = dn_ppsur
650	grid_zgr_nam%d_ppa0 = dn_ppa0
651	grid_zgr_nam%d_ppa1 = dn_ppa1
652	grid_zgr_nam%d_ppkth = dn_ppkth
653	grid_zgr_nam%d_ppacr = dn_ppacr
654	grid_zgr_nam%d_ppdzmin = dn_ppdzmin
655	grid_zgr_nam%d_pphmax = dn_pphmax
656
657	grid_zgr_nam%l_dbletanh = ln_dbletanh
658	grid_zgr_nam%d_ppa2 = dn_ppa2
659	grid_zgr_nam%d_ppkth2 = dn_ppkth2
660	grid_zgr_nam%d_ppacr2 = dn_ppacr2
661
662	grid_zgr_nam%d_e3zps_min= dn_e3zps_min
663	grid_zgr_nam%d_e3zps_rat= dn_e3zps_rat
664	! grid_zgr_nam%i_msh = in_msh
665
666	grid_zgr_nam%l_s_sh94 = ln_s_sh94
667	grid_zgr_nam%l_s_sf12 = ln_s_sf12
668	grid_zgr_nam%d_sbot_min = dn_sbot_min
669	grid_zgr_nam%d_sbot_max = dn_sbot_max
670	grid_zgr_nam%d_rmax = dn_rmax
671	grid_zgr_nam%d_hc = dn_hc
672	!
673	grid_zgr_nam%d_theta = dn_theta
674	grid_zgr_nam%d_thetb = dn_thetb
675	grid_zgr_nam%d_bb = dn_bb
676	!
677	grid_zgr_nam%l_sigcrit = ln_sigcrit
678	grid_zgr_nam%d_alpha = dn_alpha
679	grid_zgr_nam%d_efold = dn_efold
680	grid_zgr_nam%d_zs = dn_zs
681	grid_zgr_nam%d_zb_a = dn_zb_a
682	grid_zgr_nam%d_zb_b = dn_zb_b
683
684	! grid_zgr_nam%i_cla = in_cla
685
686	grid_zgr_nam%d_wdmin1 = dn_wdmin1
687	grid_zgr_nam%d_wdmin2 = dn_wdmin2
688	grid_zgr_nam%d_wdld = dn_wdld
689
690	! grid_zgr_nam%c_cfg = TRIM(cn_cfg)
691	! grid_zgr_nam%i_cfg = in_cfg
692	! grid_zgr_nam%i_bench = in_bench
693	! grid_zgr_nam%l_zoom = ln_zoom
694	grid_zgr_nam%l_c1d = ln_c1d
695	grid_zgr_nam%l_e3_dep = ln_e3_dep
696
697	! grid_zgr_nam%c_cfz = cn_cfz
698	! grid_zgr_nam%i_izoom = in_izoom
699	! grid_zgr_nam%i_jzoom = in_jzoom
700	! grid_zgr_nam%l_zoom_s = ln_zoom_s
701	! grid_zgr_nam%l_zoom_e = ln_zoom_e
702	! grid_zgr_nam%l_zoom_w = ln_zoom_w
703	! grid_zgr_nam%l_zoom_n = ln_zoom_n
704
705	ELSE
706
707	CALL logger_fatal(" GRID ZGR NAM: can't find "//TRIM(cd_namelist))
708
709	ENDIF
710
711	END FUNCTION grid_zgr_nam
712	!-------------------------------------------------------------------
713	!> @brief This subroutine fill vertical mesh
714	!>
715	!> @author J.Paul
716	!> @date September, 2015 - Initial version
717	!> @date October, 2016
718	!> - ice shelf cavity
719	!>
720	!> @param[in] td_nam
721	!> @param[in] jpi
722	!> @param[in] jpj
723	!> @param[in] jpk
724	!> @param[in] td_bathy
725	!> @param[in] td_risfdep
726	!-------------------------------------------------------------------
727	SUBROUTINE grid_zgr_fill( td_nam,jpi,jpj,jpk,td_bathy,td_risfdep )
728	IMPLICIT NONE
729	! Argument
730	TYPE(TNAMZ), INTENT(IN ) :: td_nam
731	INTEGER(i4), INTENT(IN ) :: jpi
732	INTEGER(i4), INTENT(IN ) :: jpj
733	INTEGER(i4), INTENT(IN ) :: jpk
734	TYPE(TVAR) , INTENT(INOUT) :: td_bathy
735	TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
736
737	! local variable
738	INTEGER(i4) :: il_count
739
740	REAL(dp) :: dl_bat
741
742	! loop indices
743	!----------------------------------------------------------------
744
745	CALL logger_info('GRID ZGR : vertical coordinate')
746	CALL logger_info('~~~~~~~')
747	CALL logger_info(' Namelist namzgr : set vertical coordinate')
748	CALL logger_info(' z-coordinate - full steps ln_zco = '//TRIM(fct_str(td_nam%l_zco)))
749	CALL logger_info(' z-coordinate - partial steps ln_zps = '//TRIM(fct_str(td_nam%l_zps)))
750	CALL logger_info(' s- or hybrid z-s-coordinate ln_sco = '//TRIM(fct_str(td_nam%l_sco)))
751	CALL logger_info(' ice shelf cavities ln_isfcav = '//TRIM(fct_str(td_nam%l_isfcav)))
752	CALL logger_info(' vertical scale factors ln_e3_dep = '//TRIM(fct_str(td_nam%l_e3_dep)))
753
754	il_count=0
755	IF( td_nam%l_zco ) il_count = il_count + 1
756	IF( td_nam%l_zps ) il_count = il_count + 1
757	IF( td_nam%l_sco ) il_count = il_count + 1
758	IF( il_count /= 1 )THEN
759	CALL logger_fatal(' GRID ZGR : none or several vertical coordinate options used' )
760	ENDIF
761	!
762	il_count=0
763	IF ( td_nam%l_zco .AND. td_nam%l_isfcav ) il_count = il_count + 1
764	IF ( td_nam%l_sco .AND. td_nam%l_isfcav ) il_count = il_count + 1
765	IF( il_count > 0 )THEN
766	CALL logger_fatal(' GRID ZGR : Cavity not tested/compatible with full step (zco) and sigma (ln_sco)' )
767	ENDIF
768
769	IF(.NOT. td_nam%l_e3_dep )THEN
770	CALL logger_info("Obsolescent definition of e3 scale factors is used")
771	ENDIF
772	! Build the vertical coordinate system
773	! ------------------------------------
774
775	! Reference z-coordinate system (always called)
776	CALL grid_zgr__z( td_nam,jpk )
777
778	! Bathymetry fields (levels and meters)
779	CALL grid_zgr__bat( td_nam,td_bathy,td_risfdep ) !jpi,jpj,td_bathy,td_risfdep )
780
781	! 1D config.: same bathy value over the 3x3 domain
782	IF( td_nam%l_c1d ) CALL lbc_lnk( td_bathy%d_value(:,:,1,1),'T',td_nam%i_perio,1._dp )
783
784	! z-coordinate
785	IF( td_nam%l_zco ) CALL grid_zgr__zco(jpk)
786
787	! Partial step z-coordinate
788	IF( td_nam%l_zps ) CALL grid_zgr__zps_fill( td_nam,jpi,jpj,jpk,td_bathy,td_risfdep )
789
790	! s-coordinate or hybrid z-s coordinate
791	IF( td_nam%l_sco ) CALL grid_zgr__sco_fill( td_nam,jpi,jpj,jpk,td_bathy )
792
793	! final adjustment of mbathy & check
794	! ----------------------------------
795
796	! ! correct mbathy in case of zoom subdomain
797	! IF( td_nam%l_zoom ) CALL grid_zgr__bat_zoom( td_nam,jpi,jpj )
798
799	! check bathymetry (mbathy) and suppress isolated ocean points
800	IF( .NOT. td_nam%l_c1d ) CALL grid_zgr__bat_ctl( td_nam,jpi,jpj,jpk )
801
802	! deepest ocean level for t-, u- and v-points
803	CALL grid_zgr__bot_level( ) !td_nam,jpi,jpj )
804
805	! shallowest ocean level for T-, U-, V- points
806	CALL grid_zgr__top_level( ) !td_nam,jpi,jpj )
807
808	! 1D config.: same mbathy value over the 3x3 domain
809	IF( td_nam%l_c1d ) THEN
810	dl_bat = tg_mbathy%d_value(2,2,1,1)
811	tg_mbathy%d_value(:,:,1,1) = dl_bat
812	END IF
813
814	CALL logger_info(' MIN val mbathy '//TRIM(fct_str(MINVAL( tg_mbathy%d_value(:,:,1,1) )))//&
815	& ' MAX '//TRIM(fct_str(MAXVAL( tg_mbathy%d_value(:,:,1,1) ))) )
816	CALL logger_info(' MIN val depth t '//TRIM(fct_str(MINVAL( tg_gdept_0%d_value(:,:,:,1) )))//&
817	& ' w '//TRIM(fct_str(MINVAL( tg_gdepw_0%d_value(:,:,:,1) )))//&
818	!& ' 3w '//TRIM(fct_str(MINVAL( tg_gdep3w_0%d_value(:,:,:,1) )))//&
819	& ' t '//TRIM(fct_str(MINVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
820	& ' f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
821	& ' u '//TRIM(fct_str(MINVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
822	& ' v '//TRIM(fct_str(MINVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
823	& ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
824	& ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
825	& ' w '//TRIM(fct_str(MINVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
826	CALL logger_info(' MAX val depth t '//TRIM(fct_str(MAXVAL( tg_gdept_0%d_value(:,:,:,1) )))//&
827	& ' w '//TRIM(fct_str(MAXVAL( tg_gdepw_0%d_value(:,:,:,1) ))) )!//&
828	!& ' 3w '//TRIM(fct_str(MAXVAL( tg_gdep3w_0%d_value(:,:,:,1) ))) )
829	CALL logger_info(' MAX val e3 t '//TRIM(fct_str(MAXVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
830	& ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
831	& ' u '//TRIM(fct_str(MAXVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
832	& ' v '//TRIM(fct_str(MAXVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
833	& ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
834	& ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
835	& ' w '//TRIM(fct_str(MAXVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
836
837	END SUBROUTINE grid_zgr_fill
838	!-------------------------------------------------------------------
839	!> @brief This subroutine set the depth of model levels and the resulting
840	!> vertical scale factors.
841	!>
842	!> @details
843	!>
844	!> ** Method : z-coordinate system (use in all type of coordinate)
845	!> The depth of model levels is defined from an analytical
846	!> function the derivative of which gives the scale factors.
847	!> both depth and scale factors only depend on k (1d arrays).<br/>
848	!> w-level:
849	!> - gdepw_1d = gdep(k)<br/>
850	!> - e3w_1d(k) = dk(gdep)(k) = e3(k)<br/>
851	!> t-level:
852	!> - gdept_1d = gdep(k+0.5)<br/>
853	!> - e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)<br/>
854	!>
855	!>
856	!> ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
857	!> - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
858	!>
859	!> !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
860	!>
861	!> @author J.Paul
862	!> @date September, 2015 - rewrite from zgr_z
863	!>
864	!> @param[in] td_nam
865	!> @param[in] jpk
866	!-------------------------------------------------------------------
867	SUBROUTINE grid_zgr__z(td_nam,jpk)
868	IMPLICIT NONE
869	! Argument
870	TYPE(TNAMZ), INTENT(IN ) :: td_nam
871	INTEGER(i4), INTENT(IN ) :: jpk
872
873	! local variable
874	CHARACTER(LEN=lc) :: cl_tmp
875
876	REAL(dp) :: zsur, za0, za1, zkth ! Values set from parameters in
877	REAL(dp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
878	! REAL(dp) :: zrefdep ! depth of the reference level (~10m)
879	REAL(dp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
880	REAL(dp) :: zt, zw ! temporary scalars
881	REAL(dp), PARAMETER :: dp_pp_to_be_computed = NF90_FILL_DOUBLE
882
883	! loop indices
884	INTEGER(i4) :: jk
885	!----------------------------------------------------------------
886
887	! Set variables from parameters
888	! ------------------------------
889	zkth = td_nam%d_ppkth
890	zacr = td_nam%d_ppacr
891	zdzmin = td_nam%d_ppdzmin
892	zhmax = td_nam%d_pphmax
893	zkth2 = td_nam%d_ppkth2
894	zacr2 = td_nam%d_ppacr2
895
896	! If ppa1 and ppa0 and ppsur are set to pp_to_be_computed
897	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
898	IF( td_nam%d_ppa1 == dp_pp_to_be_computed .AND. &
899	& td_nam%d_ppa0 == dp_pp_to_be_computed .AND. &
900	& td_nam%d_ppsur == dp_pp_to_be_computed ) THEN
901	!
902	za1 = ( zdzmin - zhmax / REAL(jpk-1,dp) ) &
903	& / ( TANH((1-zkth)/zacr) - zacr/REAL(jpk-1,dp) * ( LOG( COSH( (jpk - zkth) / zacr) ) &
904	& - LOG( COSH( ( 1 - zkth) / zacr) ) ) )
905	za0 = zdzmin - za1 * TANH( (1-zkth) / zacr )
906	zsur = - za0 - za1 * zacr * LOG( COSH( (1-zkth) / zacr ) )
907	! za2 ???
908	ELSE
909	za1 = td_nam%d_ppa1
910	za0 = td_nam%d_ppa0
911	zsur = td_nam%d_ppsur
912	za2 = td_nam%d_ppa2 ! optional (ldbletanh=T) double tanh parameter
913	ENDIF
914
915	CALL logger_info(' GRID ZGR Z : Reference vertical z-coordinates')
916	CALL logger_info('~~~~~~~~~~~')
917	IF( zkth == 0._dp ) THEN
918	CALL logger_info('Uniform grid with '//TRIM(fct_str(jpk-1))//' layers')
919	CALL logger_info('Total depth :'//TRIM(fct_str(zhmax)))
920	CALL logger_info('Layer thickness:'//TRIM(fct_str(zhmax/(jpk-1))))
921	ELSE
922	IF( za1 == 0._dp .AND. za0 == 0._dp .AND. zsur == 0._dp ) THEN
923	CALL logger_info('zsur, za0, za1 computed from ')
924	CALL logger_info(' zdzmin = '//TRIM(fct_str(zdzmin)))
925	CALL logger_info(' zhmax = '//TRIM(fct_str(zhmax)))
926	ENDIF
927	CALL logger_info('Value of coefficients for vertical mesh:')
928	CALL logger_info(' zsur = '//TRIM(fct_str(zsur)))
929	CALL logger_info(' za0 = '//TRIM(fct_str(za0)))
930	CALL logger_info(' za1 = '//TRIM(fct_str(za1)))
931	CALL logger_info(' zkth = '//TRIM(fct_str(zkth)))
932	CALL logger_info(' zacr = '//TRIM(fct_str(zacr)))
933	IF( td_nam%l_dbletanh ) THEN
934	CALL logger_info(' (Double tanh za2 = '//TRIM(fct_str(za2)))
935	CALL logger_info(' parameters) zkth2= '//TRIM(fct_str(zkth2)))
936	CALL logger_info(' zacr2= '//TRIM(fct_str(zacr2)))
937	ENDIF
938	ENDIF
939
940	! Reference z-coordinate (depth - scale factor at T- and W-points)
941	! ======================
942	! init
943	IF( zkth == 0._dp ) THEN ! uniform vertical grid
944	za1 = zhmax / REAL(jpk-1,dp)
945	DO jk = 1, jpk
946	zw = REAL( jk, dp )
947	zt = REAL( jk, dp ) + 0.5_dp
948	tg_gdepw_1d%d_value(1,1,jk,1) = ( zw - 1 ) * za1
949	tg_gdept_1d%d_value(1,1,jk,1) = ( zt - 1 ) * za1
950	tg_e3w_1d%d_value (1,1,jk,1) = za1
951	tg_e3t_1d%d_value (1,1,jk,1) = za1
952	END DO
953	ELSE ! Madec & Imbard 1996 function
954	IF( .NOT. td_nam%l_dbletanh ) THEN
955	DO jk = 1, jpk
956	zw = REAL( jk , dp )
957	zt = REAL( jk , dp ) + 0.5_dp
958	tg_gdepw_1d%d_value(1,1,jk,1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
959	tg_gdept_1d%d_value(1,1,jk,1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
960	tg_e3w_1d%d_value (1,1,jk,1) = za0 + za1 * TANH( (zw-zkth) / zacr )
961	tg_e3t_1d%d_value (1,1,jk,1) = za0 + za1 * TANH( (zt-zkth) / zacr )
962	END DO
963	ELSE
964	DO jk = 1, jpk
965	zw = REAL( jk, dp )
966	zt = REAL( jk, dp ) + 0.5_dp
967	! Double tanh function
968	tg_gdepw_1d%d_value(1,1,jk,1) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
969	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
970	tg_gdept_1d%d_value(1,1,jk,1) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
971	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
972	tg_e3w_1d %d_value(1,1,jk,1) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
973	& + za2 * TANH( (zw-zkth2) / zacr2 )
974	tg_e3t_1d %d_value(1,1,jk,1) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
975	& + za2 * TANH( (zt-zkth2) / zacr2 )
976	END DO
977	ENDIF
978	tg_gdepw_1d%d_value(1,1,1,1) = 0._dp ! force first w-level to be exactly at zero
979	ENDIF
980
981	IF ( td_nam%l_isfcav .OR. td_nam%l_e3_dep ) THEN
982	! need to be like this to compute the pressure gradient with ISF.
983	! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
984	! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
985	DO jk = 1, jpk-1
986	tg_e3t_1d%d_value(1,1,jk,1) = tg_gdepw_1d%d_value(1,1,jk+1,1)-tg_gdepw_1d%d_value(1,1,jk,1)
987	END DO
988	! we don't care because this level is masked in NEMO
989	tg_e3t_1d%d_value(1,1,jpk,1) = tg_e3t_1d%d_value(1,1,jpk-1,1)
990
991	DO jk = 2, jpk
992	tg_e3w_1d%d_value(1,1,jk,1) = tg_gdept_1d%d_value(1,1,jk,1) - tg_gdept_1d%d_value(1,1,jk-1,1)
993	END DO
994	tg_e3w_1d%d_value(1,1,1,1) = 2._dp * (tg_gdept_1d%d_value(1,1,1,1) - tg_gdepw_1d%d_value(1,1,1,1))
995	END IF
996
997	! unused ?
998	!!!!gm BUG in s-coordinate this does not work!
999	! ! deepest/shallowest W level Above/Below ~10m
1000	!
1001	! ! ref. depth with tolerance (10% of minimum layer thickness)
1002	! zrefdep = 10._dp - 0.1_dp * MINVAL( tg_e3w_1d%d_value(1,1,:,1) )
1003	!
1004	! ! shallowest W level Below ~10m
1005	! nlb10 = MINLOC( tg_gdepw_1d%d_value(1,1,:,1), mask = tg_gdepw_1d%d_value(1,1,:,1) > zrefdep, dim = 1 )
1006	!
1007	! ! deepest W level Above ~10m
1008	! nla10 = nlb10 - 1
1009	!!!!gm end bug
1010
1011	! control print
1012	CALL logger_info(' GRID ZGR Z : Reference z-coordinate depth and scale factors:')
1013	CALL logger_info('~~~~~~~~~~~')
1014	WRITE(cl_tmp, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
1015	CALL logger_info(cl_tmp)
1016	DO jk=1,jpk
1017	WRITE(cl_tmp, "(10x, i4, 4f9.2)" ) jk, tg_gdept_1d%d_value(1,1,jk,1), tg_gdepw_1d%d_value(1,1,jk,1), &
1018	& tg_e3t_1d%d_value (1,1,jk,1), tg_e3w_1d%d_value (1,1,jk,1)
1019	CALL logger_info(cl_tmp)
1020	ENDDO
1021
1022	! control positivity
1023	DO jk = 1, jpk
1024	IF( tg_e3w_1d%d_value (1,1,jk,1) <= 0._dp .OR. tg_e3t_1d%d_value (1,1,jk,1) <= 0._dp )THEN
1025	CALL logger_fatal( 'GRID ZGR Z: e3w_1d or e3t_1d =< 0 ' )
1026	ENDIF
1027	IF( tg_gdepw_1d%d_value(1,1,jk,1) < 0._dp .OR. tg_gdept_1d%d_value(1,1,jk,1) < 0._dp )THEN
1028	CALL logger_fatal( 'GRID ZGR Z: gdepw_1d or gdept_1d < 0 ' )
1029	ENDIF
1030	END DO
1031
1032	END SUBROUTINE grid_zgr__z
1033	!-------------------------------------------------------------------
1034	!> @brief This subroutine set bathymetry both in levels and meters
1035	!>
1036	!> @details
1037	!>
1038	!> ** Method : read or define mbathy and bathy arrays
1039	!> * level bathymetry:
1040	!> The ocean basin geometry is given by a two-dimensional array,
1041	!> mbathy, which is defined as follow :
1042	!> mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
1043	!> at t-point (ji,jj).
1044	!> = 0 over the continental t-point.
1045	!> The array mbathy is checked to verified its consistency with
1046	!> model option. in particular:
1047	!> mbathy must have at least 1 land grid-points (mbathy<=0)
1048	!> along closed boundary.
1049	!> mbathy must be cyclic IF jperio=1.
1050	!> mbathy must be lower or equal to jpk-1.
1051	!> isolated ocean grid points are suppressed from mbathy
1052	!> since they are only connected to remaining
1053	!> ocean through vertical diffusion.
1054	!> ntopo=-1 : rectangular channel or bassin with a bump
1055	!> ntopo= 0 : flat rectangular channel or basin
1056	!> ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
1057	!> bathy is read in 'bathy_meter.nc' NetCDF file
1058	!>
1059	!> ** Action : - mbathy: level bathymetry (in level index)
1060	!> - bathy : meter bathymetry (in meters)
1061	!>
1062	!> @warning do not manage case ntopo=-1 or 0
1063	!>
1064	!> @author J.Paul
1065	!> @date September, 2015 - rewrite from zgr_bat
1066	!> @date October, 2016
1067	!> - do not use anymore special case for ORCA grid.
1068	!>
1069	!> @param[in] td_nam
1070	! @param[in] jpi
1071	! @param[in] jpj
1072	!> @param[in] td_bathy
1073	!> @param[in] td_risfdep
1074	!-------------------------------------------------------------------
1075	SUBROUTINE grid_zgr__bat( td_nam,td_bathy,td_risfdep ) !jpi,jpj,td_bathy,td_risfdep )
1076	IMPLICIT NONE
1077	! Argument
1078	TYPE(TNAMZ), INTENT(IN ) :: td_nam
1079	! INTEGER(i4), INTENT(IN ) :: jpi
1080	! INTEGER(i4), INTENT(IN ) :: jpj
1081	TYPE(TVAR) , INTENT(INOUT) :: td_bathy
1082	TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
1083
1084	! local variable
1085	! INTEGER(i4) :: ii0, ii1
1086	! INTEGER(i4) :: ij0, ij1
1087
1088	REAL(dp) :: zhmin
1089
1090	! loop indices
1091	INTEGER(i4) :: jk
1092	!----------------------------------------------------------------
1093
1094	CALL logger_info(' GRID ZGR BAT : defines level and meter bathymetry')
1095	CALL logger_info(' ~~~~~~~~~~~~~')
1096	CALL logger_info(' GRID ZGR BAT : bathymetry read in file')
1097
1098	IF( td_nam%l_zco )THEN
1099	! zco : read level bathymetry
1100
1101	! read variable in bathymetry file
1102	tg_mbathy%d_value(:,:,1,1) = INT(td_bathy%d_value(:,:,1,1),i4)
1103	tg_misfdep%d_value(:,:,1,1)=1
1104	! ! =====================
1105	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. td_nam%i_cfg == 2 ) THEN ! ORCA R2 configuration
1106	! ! ! =====================
1107	! IF( td_nam%i_cla == 0 ) THEN
1108	! ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
1109	! ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
1110	! tg_mbathy%d_value(ii0:ii1,ij0:ij1,1,1) = 15
1111	! CALL logger_info('orca_r2: Gibraltar strait open at i='//&
1112	! & TRIM(fct_str(ii0))//' j='//TRIM(fct_str(ij0)) )
1113	! !
1114	! ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
1115	! ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
1116	! tg_mbathy%d_value(ii0:ii1,ij0:ij1,1,1) = 12
1117	! CALL logger_info('orca_r2: Bab el Mandeb strait open at i='//&
1118	! & TRIM(fct_str(ii0))//' j='//TRIM(fct_str(ij0)) )
1119	! ENDIF
1120	! !
1121	! ENDIF
1122
1123	ENDIF
1124
1125	IF( td_nam%l_zps .OR. td_nam%l_sco )THEN
1126	! zps or sco : read meter bathymetry
1127
1128	tg_misfdep%d_value(:,:,:,:)=1
1129
1130	IF ( td_nam%l_isfcav ) THEN
1131	WHERE( td_bathy%d_value(:,:,1,1) <= 0._dp )
1132	td_risfdep%d_value(:,:,1,1) = 0._dp
1133	END WHERE
1134
1135	! set grounded point to 0
1136	! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
1137	WHERE ( td_bathy%d_value(:,:,1,1) <= td_risfdep%d_value(:,:,1,1) + td_nam%d_isfhmin )
1138	tg_misfdep%d_value(:,:,1,1) = 0
1139	td_risfdep%d_value(:,:,1,1) = 0._dp
1140	tg_mbathy%d_value (:,:,1,1) = 0
1141	td_bathy%d_value (:,:,1,1) = 0._dp
1142	END WHERE
1143	END IF
1144	!
1145	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. td_nam%i_cfg == 2 ) THEN ! ORCA R2 configuration
1146	! !
1147	! IF( td_nam%i_cla == 0 ) THEN
1148	! ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
1149	! ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
1150	! td_bathy%d_value(ii0:ii1,ij0:ij1,1,1) = 284._dp
1151	! CALL logger_info('orca_r2: Gibraltar strait open at i='//&
1152	! & TRIM(fct_str(ii0))//' j='//TRIM(fct_str(ij0)) )
1153	! !
1154	! ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
1155	! ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
1156	! td_bathy%d_value(ii0:ii1,ij0:ij1,1,1) = 137._dp
1157	! CALL logger_info('orca_r2: Bab el Mandeb strait open at i='//&
1158	! & TRIM(fct_str(ii0))//' j='//TRIM(fct_str(ij0)) )
1159	! ENDIF
1160	! !
1161	! ENDIF
1162	!
1163	ENDIF
1164
1165	!== NO closed seas or lakes ==!
1166	! already done
1167
1168	IF ( .NOT. td_nam%l_sco ) THEN
1169	!== set a minimum depth ==!
1170
1171	IF( td_nam%d_hmin < 0._dp ) THEN
1172	! from a nb of level
1173	jk = - INT(td_nam%d_hmin, i4)
1174	ELSE
1175	! from a depth
1176	jk = MINLOC( tg_gdepw_1d%d_value(1,1,:,1), &
1177	& MASK = tg_gdepw_1d%d_value(1,1,:,1) > td_nam%d_hmin, &
1178	& DIM = 1)
1179	ENDIF
1180
1181	! minimum depth = ik+1 w-levels
1182	zhmin = tg_gdepw_1d%d_value(1,1,jk+1,1)
1183	WHERE( td_bathy%d_value(:,:,1,1) <= 0._dp )
1184	! min=0 over the lands
1185	td_bathy%d_value(:,:,1,1) = 0._dp
1186	ELSE WHERE
1187	! min=zhmin over the oceans
1188	td_bathy%d_value(:,:,1,1) = MAX(zhmin, td_bathy%d_value(:,:,1,1))
1189	END WHERE
1190	CALL logger_info('GRID ZGR BAT: Minimum ocean depth: '//&
1191	& TRIM(fct_str(zhmin))//' minimum number of ocean'//&
1192	& ' levels : '//TRIM(fct_str(jk)))
1193	ENDIF
1194
1195	END SUBROUTINE grid_zgr__bat
1196	!-------------------------------------------------------------------
1197	!> @brief This subroutine define the z-coordinate system
1198	!>
1199	!> @details
1200	!> set 3D coord. arrays to reference 1D array
1201	!>
1202	!> @author J.Paul
1203	!> @date September, 2015 - rewrite from zgr_zco
1204	!>
1205	!> @param[in] jpk
1206	!-------------------------------------------------------------------
1207	SUBROUTINE grid_zgr__zco(jpk)
1208	IMPLICIT NONE
1209	! Argument
1210	INTEGER(i4), INTENT(IN ) :: jpk
1211
1212	! local variable
1213	! loop indices
1214	INTEGER(i4) :: jk
1215	!----------------------------------------------------------------
1216
1217	DO jk = 1, jpk
1218	tg_gdept_0%d_value (:,:,jk,1) = tg_gdept_1d%d_value(1,1,jk,1)
1219	tg_gdepw_0%d_value (:,:,jk,1) = tg_gdepw_1d%d_value(1,1,jk,1)
1220	!tg_gdep3w_0%d_value(:,:,jk,1) = tg_gdepw_1d%d_value(1,1,jk,1)
1221	tg_e3t_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value (1,1,jk,1)
1222	tg_e3u_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value (1,1,jk,1)
1223	tg_e3v_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value (1,1,jk,1)
1224	tg_e3f_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value (1,1,jk,1)
1225	tg_e3w_0%d_value (:,:,jk,1) = tg_e3w_1d%d_value (1,1,jk,1)
1226	tg_e3uw_0%d_value (:,:,jk,1) = tg_e3w_1d%d_value (1,1,jk,1)
1227	tg_e3vw_0%d_value (:,:,jk,1) = tg_e3w_1d%d_value (1,1,jk,1)
1228	END DO
1229
1230	END SUBROUTINE grid_zgr__zco
1231	! !-------------------------------------------------------------------
1232	! !> @brief This subroutine :
1233	! !> - close zoom domain boundary if necessary
1234	! !> - suppress Med Sea from ORCA R2 and R05 arctic zoom
1235	! !>
1236	! !> @author J.Paul
1237	! !> @date September, 2015 - Initial version
1238	! !>
1239	! !> @param[in] td_nam
1240	! !> @param[in] jpi
1241	! !> @param[in] jpj
1242	! !-------------------------------------------------------------------
1243	! SUBROUTINE grid_zgr__bat_zoom(td_nam,jpi,jpj)
1244	! IMPLICIT NONE
1245	! ! Argument
1246	! TYPE(TNAMZ), INTENT(IN ) :: td_nam
1247	! INTEGER(i4), INTENT(IN ) :: jpi
1248	! INTEGER(i4), INTENT(IN ) :: jpj
1249	!
1250	! ! local variable
1251	! INTEGER(i4) :: jpizoom
1252	! INTEGER(i4) :: jpjzoom
1253	!
1254	! INTEGER(i4) :: ii0, ii1
1255	! INTEGER(i4) :: ij0, ij1
1256	! ! loop indices
1257	! !----------------------------------------------------------------
1258	!
1259	! CALL logger_info('GRID ZGR BAT ZOOM : modify the level bathymetry for zoom domain')
1260	! CALL logger_info('~~~~~~~~~~~~')
1261	!
1262	! jpizoom=td_nam%i_izoom
1263	! jpjzoom=td_nam%i_jzoom
1264	!
1265	! ! Forced closed boundary if required
1266	! IF( td_nam%l_zoom_s ) tg_mbathy%d_value( : , jpjzoom ,1,1) = 0
1267	! IF( td_nam%l_zoom_w ) tg_mbathy%d_value( jpizoom , : ,1,1) = 0
1268	! IF( td_nam%l_zoom_e ) tg_mbathy%d_value( jpi+jpizoom-1, : ,1,1) = 0
1269	! IF( td_nam%l_zoom_n ) tg_mbathy%d_value( : ,jpj+jpjzoom-1,1,1) = 0
1270	!
1271	! ! Configuration specific domain modifications
1272	! ! (here, ORCA arctic configuration: suppress Med Sea)
1273	! IF( TRIM(td_nam%c_cfg) == "orca" .AND. &
1274	! & TRIM(td_nam%c_cfz) == "arctic" ) THEN
1275	! SELECT CASE ( td_nam%i_cfg )
1276	! ! ! =======================
1277	! CASE ( 2 ) ! ORCA_R2 configuration
1278	! ! ! =======================
1279	! CALL logger_info('ORCA R2 arctic zoom: suppress the Med Sea')
1280	! ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
1281	! ij0 = 98 ; ij1 = 110
1282	! ! ! =======================
1283	! CASE ( 05 ) ! ORCA_R05 configuration
1284	! ! ! =======================
1285	! CALL logger_info('ORCA R05 arctic zoom: suppress the Med Sea')
1286	! ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
1287	! ij0 = 314 ; ij1 = 370
1288	! END SELECT
1289	! !
1290	! tg_mbathy%d_value( ii0:ii1, ij0:ij1, 1, 1) = 0 ! zero over the Med Sea boxe
1291	! !
1292	! ENDIF
1293	!
1294	! END SUBROUTINE grid_zgr__bat_zoom
1295	!-------------------------------------------------------------------
1296	!> @brief This subroutine check the bathymetry in levels
1297	!>
1298	!> @details
1299	!>
1300	!>
1301	!> ** Method : The array mbathy is checked to verified its consistency
1302	!> with the model options. in particular:
1303	!> mbathy must have at least 1 land grid-points (mbathy<=0)
1304	!> along closed boundary.
1305	!> mbathy must be cyclic IF jperio=1.
1306	!> mbathy must be lower or equal to jpk-1.
1307	!> isolated ocean grid points are suppressed from mbathy
1308	!> since they are only connected to remaining
1309	!> ocean through vertical diffusion.
1310	!> C A U T I O N : mbathy will be modified during the initializa-
1311	!> tion phase to become the number of non-zero w-levels of a water
1312	!> column, with a minimum value of 1.
1313	!>
1314	!> ** Action : - update mbathy: level bathymetry (in level index)
1315	!> - update bathy : meter bathymetry (in meters)
1316	!>
1317	!> @author J.Paul
1318	!> @date September, 2015 - Initial version
1319	!>
1320	!> @param[in] td_nam
1321	!> @param[in] jpi
1322	!> @param[in] jpj
1323	!> @param[in] jpk
1324	!-------------------------------------------------------------------
1325	SUBROUTINE grid_zgr__bat_ctl(td_nam,jpi,jpj,jpk)
1326	IMPLICIT NONE
1327	! Argument
1328	TYPE(TNAMZ), INTENT(IN ) :: td_nam
1329	INTEGER(i4), INTENT(IN ) :: jpi
1330	INTEGER(i4), INTENT(IN ) :: jpj
1331	INTEGER(i4), INTENT(IN ) :: jpk
1332
1333	! local variable
1334	INTEGER(i4) :: icompt, ibtest, ikmax ! temporary integers
1335
1336	! loop indices
1337	INTEGER(i4) :: ji
1338	INTEGER(i4) :: jj
1339	INTEGER(i4) :: jl
1340	!----------------------------------------------------------------
1341
1342	CALL logger_info('GRID ZGR BAT CTL: check the bathymetry')
1343	CALL logger_info('~~~~~~~~~~~~~~~')
1344	CALL logger_info(' suppress isolated ocean grid points')
1345	CALL logger_info(' -----------------------------------')
1346
1347	icompt = 0
1348	DO jl = 1, 2
1349	IF( td_nam%i_perio == 1 .OR. &
1350	& td_nam%i_perio == 4 .OR. &
1351	& td_nam%i_perio == 6 ) THEN
1352	tg_mbathy%d_value( 1 ,:,1,1) = tg_mbathy%d_value(jpi-1,:,1,1) ! local domain is cyclic east-west
1353	tg_mbathy%d_value(jpi,:,1,1) = tg_mbathy%d_value( 2 ,:,1,1)
1354	ENDIF
1355	DO jj = 2, jpj-1
1356	DO ji = 2, jpi-1
1357	ibtest = MAX( tg_mbathy%d_value(ji-1,jj ,1,1), &
1358	& tg_mbathy%d_value(ji+1,jj ,1,1), &
1359	& tg_mbathy%d_value(ji ,jj-1,1,1), &
1360	& tg_mbathy%d_value(ji ,jj+1,1,1) )
1361	IF( ibtest < tg_mbathy%d_value(ji,jj,1,1) ) THEN
1362	CALL logger_info(' the number of ocean level at '//&
1363	& 'grid-point (i,j) = ('//TRIM(fct_str(ji))//&
1364	& ','//TRIM(fct_str(jj))//') is changed from '//&
1365	& TRIM(fct_str(tg_mbathy%d_value(ji,jj,1,1)))//' to '//&
1366	& TRIM(fct_str(ibtest)) )
1367	tg_mbathy%d_value(ji,jj,1,1) = ibtest
1368	icompt = icompt + 1
1369	ENDIF
1370	END DO
1371	END DO
1372	END DO
1373
1374	!! lk_mpp not used
1375
1376	IF( icompt == 0 ) THEN
1377	CALL logger_info(' no isolated ocean grid points')
1378	ELSE
1379	CALL logger_info(TRIM(fct_str(icompt))//' ocean grid points suppressed')
1380	ENDIF
1381
1382	!! lk_mpp not used
1383
1384	! ! East-west cyclic boundary conditions
1385	IF( td_nam%i_perio == 0 ) THEN
1386	CALL logger_info(' mbathy set to 0 along east and west boundary:'//&
1387	& ' nperio = '//TRIM(fct_str(td_nam%i_perio)) )
1388	!! lk_mpp not used
1389	IF( td_nam%l_zco .OR. td_nam%l_zps ) THEN
1390	tg_mbathy%d_value( 1 ,:,1,1) = 0
1391	tg_mbathy%d_value(jpi,:,1,1) = 0
1392	ELSE
1393	tg_mbathy%d_value( 1 ,:,1,1) = jpk-1
1394	tg_mbathy%d_value(jpi,:,1,1) = jpk-1
1395	ENDIF
1396	ELSEIF( td_nam%i_perio == 1 .OR. &
1397	& td_nam%i_perio == 4 .OR. &
1398	& td_nam%i_perio == 6 ) THEN
1399	CALL logger_info(' east-west cyclic boundary conditions on mbathy:'//&
1400	& ' nperio = '//TRIM(fct_str(td_nam%i_perio)) )
1401	tg_mbathy%d_value( 1 ,:,1,1) = tg_mbathy%d_value(jpi-1,:,1,1)
1402	tg_mbathy%d_value(jpi,:,1,1) = tg_mbathy%d_value( 2 ,:,1,1)
1403	ELSEIF( td_nam%i_perio == 2 ) THEN
1404	CALL logger_info(' equatorial boundary conditions on mbathy:'//&
1405	' nperio = '//TRIM(fct_str(td_nam%i_perio)) )
1406	ELSE
1407	CALL logger_info(' e r r o r')
1408	CALL logger_info(' parameter , nperio = '//TRIM(fct_str(td_nam%i_perio)) )
1409	! STOP 'dom_mba'
1410	ENDIF
1411
1412	! Boundary condition on mbathy
1413	!!gm !!bug ??? think about it !
1414	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
1415	CALL lbc_lnk( tg_mbathy%d_value(:,:,1,1), 'T', td_nam%i_perio, 1._dp )
1416
1417
1418	! Number of ocean level inferior or equal to jpkm1
1419	ikmax = 0
1420	DO jj = 1, jpj
1421	DO ji = 1, jpi
1422	ikmax = MAX( ikmax, INT(tg_mbathy%d_value(ji,jj,1,1),i4) )
1423	END DO
1424	END DO
1425	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
1426	IF( ikmax > jpk-1 ) THEN
1427	CALL logger_info(' maximum number of ocean level = '//TRIM(fct_str(ikmax))//' > jpk-1')
1428	CALL logger_info(' change jpk to '//TRIM(fct_str(ikmax+1))//' to use the exact ead bathymetry')
1429	ELSE IF( ikmax < jpk-1 ) THEN
1430	CALL logger_info(' maximum number of ocean level = '//TRIM(fct_str(ikmax))//' < jpk-1')
1431	CALL logger_info(' you can decrease jpk to '//TRIM(fct_str(ikmax+1)))
1432	ENDIF
1433
1434	END SUBROUTINE grid_zgr__bat_ctl
1435	!-------------------------------------------------------------------
1436	!> @brief This subroutine defines the vertical index of ocean bottom (mbk. arrays)
1437	!>
1438	!> @details
1439	!>
1440	!> ** Method : computes from mbathy with a minimum value of 1 over land
1441	!>
1442	!> ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
1443	!> ocean level at t-, u- & v-points
1444	!> (min value = 1 over land)
1445	!>
1446	!>
1447	!> @author J.Paul
1448	!> @date September, 2015 - rewrite from zgr_bot_level
1449	!>
1450	! @param[in] td_nam
1451	! @param[in] jpi
1452	! @param[in] jpj
1453	!-------------------------------------------------------------------
1454	SUBROUTINE grid_zgr__bot_level( )!td_nam,jpi,jpj)
1455	IMPLICIT NONE
1456	! Argument
1457	! TYPE(TNAMZ), INTENT(IN ) :: td_nam
1458	! INTEGER(i4), INTENT(IN ) :: jpi
1459	! INTEGER(i4), INTENT(IN ) :: jpj
1460
1461	! local variable
1462
1463	! loop indices
1464	! INTEGER(i4) :: ji
1465	! INTEGER(i4) :: jj
1466	!----------------------------------------------------------------
1467
1468	CALL logger_info('GRID ZGR BOT LEVEL: ocean bottom k-index of T-, U-, V- and W-levels ')
1469	CALL logger_info(' ~~~~~~~~~~~~~')
1470
1471	! bottom k-index of T-level (=1 over land)
1472	tg_mbkt%d_value(:,:,1,1) = MAX( tg_mbathy%d_value(:,:,1,1) , 1._dp )
1473
1474	! ! bottom k-index of W-level = mbkt+1
1475	! ! bottom k-index of u- (v-) level
1476	! DO jj = 1, jpj-1 ! bottom k-index of u- (v-) level
1477	! DO ji = 1, jpi-1
1478	! tg_mbku%d_value(ji,jj,1,1) = MIN( tg_mbkt%d_value(ji+1,jj ,1,1) , &
1479	! & tg_mbkt%d_value(ji ,jj ,1,1) )
1480	!
1481	! tg_mbkv%d_value(ji,jj,1,1) = MIN( tg_mbkt%d_value(ji ,jj+1,1,1) , &
1482	! & tg_mbkt%d_value(ji ,jj ,1,1) )
1483	! END DO
1484	! END DO
1485	!
1486	! CALL lbc_lnk(tg_mbku%d_value(:,:,1,1),'U', td_nam%i_perio, 1._dp)
1487	! CALL lbc_lnk(tg_mbkv%d_value(:,:,1,1),'U', td_nam%i_perio, 1._dp)
1488
1489	END SUBROUTINE grid_zgr__bot_level
1490	!-------------------------------------------------------------------
1491	!> @brief This subroutine defines the vertical index of ocean top (mik. arrays)
1492	!>
1493	!> @details
1494	!>
1495	!> ** Method : computes from misfdep with a minimum value of 1
1496	!>
1497	!> ** Action : mikt, miku, mikv : vertical indices of the shallowest
1498	!> ocean level at t-, u- & v-points
1499	!> (min value = 1)
1500	!>
1501	!> @author J.Paul
1502	!> @date September, 2015 - rewrite from zgr_top_level
1503	!>
1504	! @param[in] td_nam
1505	! @param[in] jpi
1506	! @param[in] jpj
1507	!-------------------------------------------------------------------
1508	SUBROUTINE grid_zgr__top_level( )!td_nam,jpi,jpj)
1509	IMPLICIT NONE
1510	! Argument
1511	! TYPE(TNAMZ), INTENT(IN ) :: td_nam
1512	! INTEGER(i4), INTENT(IN ) :: jpi
1513	! INTEGER(i4), INTENT(IN ) :: jpj
1514	! local variable
1515
1516	! loop indices
1517	! INTEGER(i4) :: ji
1518	! INTEGER(i4) :: jj
1519	!----------------------------------------------------------------
1520
1521	CALL logger_info(' GRID ZGR TOP LEVEL : ocean top k-index of T-, U-, V- and W-levels ')
1522	CALL logger_info(' ~~~~~~~~~~~~~')
1523
1524	! top k-index of T-level (=1)
1525	tg_mikt%d_value(:,:,1,1) = MAX( tg_misfdep%d_value(:,:,1,1) , 1._dp )
1526
1527	! ! top k-index of W-level (=mikt)
1528	! ! top k-index of U- (U-) level
1529	! DO jj = 1, jpj-1 ! top k-index of U- (U-) level
1530	! DO ji = 1, jpi-1
1531	! tg_miku%d_value(ji,jj,1,1) = MAX( tg_mikt%d_value(ji+1,jj ,1,1), &
1532	! & tg_mikt%d_value(ji ,jj ,1,1) )
1533	!
1534	! tg_mikv%d_value(ji,jj,1,1) = MAX( tg_mikt%d_value(ji ,jj+1,1,1), &
1535	! & tg_mikt%d_value(ji ,jj ,1,1) )
1536	!
1537	! tg_mikf%d_value(ji,jj,1,1) = MAX( tg_mikt%d_value(ji ,jj+1,1,1), &
1538	! & tg_mikt%d_value(ji ,jj ,1,1), &
1539	! & tg_mikt%d_value(ji+1,jj ,1,1), &
1540	! & tg_mikt%d_value(ji+1,jj+1,1,1) )
1541	! END DO
1542	! END DO
1543	!
1544	! CALL lbc_lnk(tg_miku%d_value(:,:,1,1),'U',td_nam%i_perio,1._dp)
1545	! CALL lbc_lnk(tg_mikv%d_value(:,:,1,1),'V',td_nam%i_perio,1._dp)
1546	! CALL lbc_lnk(tg_mikf%d_value(:,:,1,1),'F',td_nam%i_perio,1._dp)
1547
1548	END SUBROUTINE grid_zgr__top_level
1549	!-------------------------------------------------------------------
1550	!> @brief This function initialise global variable needed to compute vertical
1551	!> mesh
1552	!>
1553	!> @author J.Paul
1554	!> @date September, 2015 - Initial version
1555	!>
1556	!> @param[in] jpi
1557	!> @param[in] jpj
1558	!-------------------------------------------------------------------
1559	SUBROUTINE grid_zgr_zps_init(jpi,jpj)
1560	IMPLICIT NONE
1561	! Argument
1562	INTEGER(i4), INTENT(IN) :: jpi
1563	INTEGER(i4), INTENT(IN) :: jpj
1564
1565	! local variable
1566	REAL(dp), DIMENSION(jpi,jpj) :: dl_tmp
1567
1568	! loop indices
1569	!----------------------------------------------------------------
1570
1571	dl_tmp(:,:)=dp_fill
1572
1573	tg_e3tp =var_init('e3t_ps ',dl_tmp(:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
1574	tg_e3wp =var_init('e3w_ps ',dl_tmp(:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
1575
1576	END SUBROUTINE grid_zgr_zps_init
1577	!-------------------------------------------------------------------
1578	!> @brief This function clean hgr structure
1579	!>
1580	!> @author J.Paul
1581	!> @date September, 2015 - Initial version
1582	!>
1583	!-------------------------------------------------------------------
1584	SUBROUTINE grid_zgr_zps_clean()
1585	IMPLICIT NONE
1586	! Argument
1587
1588	! local variable
1589	! loop indices
1590	!----------------------------------------------------------------
1591
1592	CALL var_clean(tg_e3tp )
1593	CALL var_clean(tg_e3wp )
1594
1595	END SUBROUTINE grid_zgr_zps_clean
1596	!-------------------------------------------------------------------
1597	!> @brief This subroutine define the depth and vertical scale factor in partial step
1598	!> z-coordinate case
1599	!>
1600	!> @details
1601	!> ** Method : Partial steps : computes the 3D vertical scale factors
1602	!> of T-, U-, V-, W-, UW-, VW and F-points that are associated with
1603	!> a partial step representation of bottom topography.
1604	!>
1605	!> The reference depth of model levels is defined from an analytical
1606	!> function the derivative of which gives the reference vertical
1607	!> scale factors.
1608	!> From depth and scale factors reference, we compute there new value
1609	!> with partial steps on 3d arrays ( i, j, k ).
1610	!>
1611	!> w-level: gdepw_0(i,j,k) = gdep(k)
1612	!> e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
1613	!> t-level: gdept_0(i,j,k) = gdep(k+0.5)
1614	!> e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
1615	!>
1616	!> With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
1617	!> we find the mbathy index of the depth at each grid point.
1618	!> This leads us to three cases:
1619	!>
1620	!> - bathy = 0 => mbathy = 0
1621	!> - 1 < mbathy < jpkm1
1622	!> - bathy > gdepw_0(jpk) => mbathy = jpkm1
1623	!>
1624	!> Then, for each case, we find the new depth at t- and w- levels
1625	!> and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
1626	!> and f-points.
1627	!>
1628	!> This routine is given as an example, it must be modified
1629	!> following the user s desiderata. nevertheless, the output as
1630	!> well as the way to compute the model levels and scale factors
1631	!> must be respected in order to insure second order accuracy
1632	!> schemes.
1633	!>
1634	!> @warrning gdept_1d, gdepw_1d and e3._1d are positives
1635	!> gdept_0, gdepw_0 and e3. are positives
1636	!>
1637	!> Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
1638	!> set 3D coord. arrays to reference 1D array
1639	!>
1640	!> @author J.Paul
1641	!> @date September, 2015 - rewrite from zgr_zps
1642	!>
1643	!> @param[in] td_nam
1644	!> @param[in] jpi
1645	!> @param[in] jpj
1646	!> @param[in] jpk
1647	!> @param[inout] td_bathy
1648	!> @param[inout] td_risfdep
1649	!-------------------------------------------------------------------
1650	SUBROUTINE grid_zgr__zps_fill( td_nam, jpi,jpj,jpk,td_bathy, td_risfdep )
1651	IMPLICIT NONE
1652	! Argument
1653	TYPE(TNAMZ), INTENT(IN ) :: td_nam
1654	INTEGER(i4), INTENT(IN ) :: jpi
1655	INTEGER(i4), INTENT(IN ) :: jpj
1656	INTEGER(i4), INTENT(IN ) :: jpk
1657	TYPE(TVAR) , INTENT(INOUT) :: td_bathy
1658	TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
1659
1660	! local variable
1661	REAL(dp) :: zmax ! Maximum depth
1662	REAL(dp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
1663	REAL(dp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1664	REAL(dp) :: zdiff ! temporary scalar
1665
1666	! loop indices
1667	INTEGER(i4) :: ji
1668	INTEGER(i4) :: jj
1669	INTEGER(i4) :: jk
1670
1671	INTEGER(i4) :: ik
1672	INTEGER(i4) :: it
1673	!----------------------------------------------------------------
1674
1675	CALL logger_info('GRID ZGR ZPS : z-coordinate with partial steps')
1676	CALL logger_info('mbathy is recomputed : bathy_level file is NOT used')
1677
1678	! bathymetry in level (from bathy_meter)
1679
1680	! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1681	zmax = tg_gdepw_1d%d_value(1,1,jpk,1) + tg_e3t_1d%d_value(1,1,jpk,1)
1682
1683	! bounded value of bathy (min already set at the end of zgr_bat)
1684	td_bathy%d_value(:,:,1,1) = MIN(zmax, td_bathy%d_value(:,:,1,1))
1685
1686	WHERE( td_bathy%d_value(:,:,1,1) == 0._dp )
1687	! land : set mbathy to 0
1688	tg_mbathy%d_value(:,:,1,1) = 0
1689	ELSE WHERE
1690	! ocean : initialize mbathy to the max ocean level
1691	tg_mbathy%d_value(:,:,1,1) = jpk-1
1692	END WHERE
1693
1694	! Compute mbathy for ocean points (i.e. the number of ocean levels)
1695	! find the number of ocean levels such that the last level thickness
1696	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1697	! e3t_1d is the reference level thickness
1698	DO jk = jpk-1, 1, -1
1699	zdepth = tg_gdepw_1d%d_value(1,1,jk,1) + MIN( td_nam%d_e3zps_min, &
1700	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,jk,1))
1701
1702	WHERE( 0._dp < td_bathy%d_value(:,:,1,1) .AND. &
1703	& td_bathy%d_value(:,:,1,1) <= zdepth )
1704	tg_mbathy%d_value(:,:,1,1) = jk-1
1705	END WHERE
1706	END DO
1707
1708	! Scale factors and depth at T- and W-points
1709	DO jk = 1, jpk
1710	! intitialization to the reference z-coordinate
1711	tg_gdept_0%d_value(:,:,jk,1) = tg_gdept_1d%d_value(1,1,jk,1)
1712	tg_gdepw_0%d_value(:,:,jk,1) = tg_gdepw_1d%d_value(1,1,jk,1)
1713	tg_e3t_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value (1,1,jk,1)
1714	tg_e3w_0%d_value (:,:,jk,1) = tg_e3w_1d%d_value (1,1,jk,1)
1715	END DO
1716
1717	IF ( td_nam%l_isfcav )THEN
1718	! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
1719	CALL grid_zgr__isf_fill( td_nam, jpi,jpj,jpk, td_bathy, td_risfdep )
1720
1721	ELSE ! .NOT. td_nam%l_isfcav
1722	!
1723	DO jj = 1, jpj
1724	DO ji = 1, jpi
1725	ik = tg_mbathy%d_value(ji,jj,1,1)
1726	IF( ik > 0 ) THEN
1727	! ocean point only
1728	IF( ik == jpk-1 ) THEN
1729	! max ocean level case
1730	zdepwp = td_bathy%d_value(ji,jj,1,1)
1731	ze3tp = td_bathy%d_value(ji,jj,1,1) - tg_gdepw_1d%d_value(1,1,ik,1)
1732	ze3wp = 0.5_dp * tg_e3w_1d%d_value(1,1,ik,1) &
1733	& * ( 1._dp + ( ze3tp/tg_e3t_1d%d_value(1,1,ik,1) ) )
1734	tg_e3t_0%d_value (ji,jj,ik ,1) = ze3tp
1735	tg_e3t_0%d_value (ji,jj,ik+1,1) = ze3tp
1736	tg_e3w_0%d_value (ji,jj,ik ,1) = ze3wp
1737	tg_e3w_0%d_value (ji,jj,ik+1,1) = ze3tp
1738
1739	tg_gdepw_0%d_value(ji,jj,ik+1,1) = zdepwp
1740	tg_gdept_0%d_value(ji,jj,ik ,1) = tg_gdept_1d%d_value( 1, 1,ik-1,1) + ze3wp
1741	tg_gdept_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value (ji,jj,ik ,1) + ze3tp
1742	!
1743	ELSE
1744	! standard case
1745	IF( td_bathy%d_value(ji,jj,1,1) <= tg_gdepw_1d%d_value(1,1,ik+1,1) ) THEN
1746	tg_gdepw_0%d_value(ji,jj,ik+1,1) = td_bathy%d_value(ji,jj,1,1)
1747	ELSE
1748	tg_gdepw_0%d_value(ji,jj,ik+1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
1749	ENDIF
1750	!gm Bug? check the gdepw_1d
1751	! ... on ik
1752	tg_gdept_0%d_value(ji,jj,ik,1) = tg_gdepw_1d%d_value( 1, 1,ik ,1) + &
1753	& ( tg_gdepw_0%d_value(ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
1754	& * ( (tg_gdept_1d%d_value( 1, 1,ik ,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
1755	& / (tg_gdepw_1d%d_value( 1, 1,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) )
1756
1757	tg_e3t_0%d_value (ji,jj,ik,1) = tg_e3t_1d%d_value ( 1, 1,ik ,1) &
1758	& * ( tg_gdepw_0%d_value (ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
1759	& / ( tg_gdepw_1d%d_value( 1, 1,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) )
1760
1761	tg_e3w_0%d_value (ji,jj,ik,1) = 0.5_dp &
1762	& * ( tg_gdepw_0%d_value(ji,jj,ik+1,1) &
1763	& + tg_gdepw_1d%d_value(1,1 ,ik+1,1) &
1764	& - tg_gdepw_1d%d_value(1,1 ,ik ,1) * 2._dp ) &
1765	& * ( tg_e3w_1d%d_value (1 ,1 ,ik ,1) &
1766	& / (tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1)) )
1767
1768	! ... on ik+1
1769	tg_e3w_0%d_value (ji,jj,ik+1,1) = tg_e3t_0%d_value (ji,jj,ik,1)
1770	tg_e3t_0%d_value (ji,jj,ik+1,1) = tg_e3t_0%d_value (ji,jj,ik,1)
1771	tg_gdept_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value(ji,jj,ik,1) + tg_e3t_0%d_value(ji,jj,ik,1)
1772	ENDIF
1773	ENDIF
1774	END DO
1775	END DO
1776	!
1777	it = 0
1778	DO jj = 1, jpj
1779	DO ji = 1, jpi
1780	ik = tg_mbathy%d_value(ji,jj,1,1)
1781	IF( ik > 0 ) THEN ! ocean point only
1782	tg_e3tp%d_value (ji,jj,1,1) = tg_e3t_0%d_value(ji,jj,ik,1)
1783	tg_e3wp%d_value (ji,jj,1,1) = tg_e3w_0%d_value(ji,jj,ik,1)
1784	! test
1785	zdiff= tg_gdepw_0%d_value(ji,jj,ik+1,1) - tg_gdept_0%d_value(ji,jj,ik,1)
1786	IF( zdiff <= 0._dp ) THEN
1787	it = it + 1
1788	CALL logger_info(' it = '//TRIM(fct_str(it))//&
1789	& ' ik = '//TRIM(fct_str(ik))//&
1790	& ' (i,j) = ('//TRIM(fct_str(ji))//','//TRIM(fct_str(jj))//')')
1791	CALL logger_info(' bathy = '//TRIM(fct_str(td_bathy%d_value(ji,jj,1,1))))
1792	CALL logger_info(' gdept_0 = '//TRIM(fct_str( tg_gdept_0%d_value(ji,jj,ik ,1) ))//&
1793	& ' gdepw_0 = '//TRIM(fct_str( tg_gdepw_0%d_value(ji,jj,ik+1,1) ))//&
1794	& ' zdiff = '//TRIM(fct_str(zdiff)) )
1795	CALL logger_info(' e3tp = '//TRIM(fct_str(tg_e3t_0%d_value(ji,jj,ik,1)))//&
1796	& ' e3wp = '//TRIM(fct_str(tg_e3w_0%d_value(ji,jj,ik,1))) )
1797	ENDIF
1798	ENDIF
1799	END DO
1800	END DO
1801	ENDIF
1802	!
1803	! IF ( td_nam%l_isfcav ) THEN
1804	! ! (ISF) Definition of e3t, u, v, w for ISF case
1805	! CALL grid_zgr__isf_fill_e3x( jpi,jpj, &
1806	! & td_risfdep )
1807	! END IF
1808
1809	! Scale factors and depth at U-, V-, UW and VW-points
1810	DO jk = 1, jpk
1811	! initialisation to z-scale factors
1812	tg_e3u_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1813	tg_e3v_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1814	tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
1815	tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
1816	END DO
1817
1818	! Computed as the minimum of neighbooring scale factors
1819	DO jk = 1,jpk
1820	DO jj = 1, jpj - 1
1821	DO ji = 1, jpi - 1
1822	tg_e3u_0%d_value (ji,jj,jk,1) = MIN( tg_e3t_0%d_value(ji,jj,jk,1), tg_e3t_0%d_value(ji+1,jj ,jk,1) )
1823	tg_e3v_0%d_value (ji,jj,jk,1) = MIN( tg_e3t_0%d_value(ji,jj,jk,1), tg_e3t_0%d_value(ji ,jj+1,jk,1) )
1824	tg_e3uw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji+1,jj ,jk,1) )
1825	tg_e3vw_0%d_value(ji,jj,jk,1) = MIN( tg_e3w_0%d_value(ji,jj,jk,1), tg_e3w_0%d_value(ji ,jj+1,jk,1) )
1826	END DO
1827	END DO
1828	END DO
1829
1830	IF ( td_nam%l_isfcav ) THEN
1831	! (ISF) define e3uw (adapted for 2 cells in the water column)
1832	CALL grid_zgr__isf_fill_e3uw(jpi,jpj)
1833	END IF
1834
1835	! lateral boundary conditions
1836	CALL lbc_lnk( tg_e3u_0%d_value (:,:,:,1), 'U', td_nam%i_perio, 1._dp )
1837	CALL lbc_lnk( tg_e3v_0%d_value (:,:,:,1), 'V', td_nam%i_perio, 1._dp )
1838	CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
1839	CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
1840
1841	! set to z-scale factor if zero (i.e. along closed boundaries)
1842	DO jk = 1, jpk
1843	WHERE( tg_e3u_0%d_value (:,:,jk,1) == 0._dp ) tg_e3u_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1844	WHERE( tg_e3v_0%d_value (:,:,jk,1) == 0._dp ) tg_e3v_0%d_value (:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1845	WHERE( tg_e3uw_0%d_value(:,:,jk,1) == 0._dp ) tg_e3uw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
1846	WHERE( tg_e3vw_0%d_value(:,:,jk,1) == 0._dp ) tg_e3vw_0%d_value(:,:,jk,1) = tg_e3w_1d%d_value(1,1,jk,1)
1847	END DO
1848
1849	!! Scale factor at F-point
1850	!DO jk = 1, jpk
1851	! ! initialisation to z-scale factors
1852	! tg_e3f_0%d_value(:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1853	!END DO
1854	!
1855	!! Computed as the minimum of neighbooring V-scale factors
1856	!DO jk = 1, jpk
1857	! DO jj = 1, jpj - 1
1858	! DO ji = 1, jpi - 1
1859	! tg_e3f_0%d_value(ji,jj,jk,1) = MIN( tg_e3v_0%d_value(ji,jj,jk,1), tg_e3v_0%d_value(ji+1,jj,jk,1) )
1860	! END DO
1861	! END DO
1862	!END DO
1863	!! Lateral boundary conditions
1864	!CALL lbc_lnk( tg_e3f_0%d_value(:,:,:,1), 'F', td_nam%i_perio, 1._dp )
1865	!
1866	!! set to z-scale factor if zero (i.e. along closed boundaries)
1867	!DO jk = 1, jpk
1868	! WHERE( tg_e3f_0%d_value(:,:,jk,1) == 0._dp ) tg_e3f_0%d_value(:,:,jk,1) = tg_e3t_1d%d_value(1,1,jk,1)
1869	!END DO
1870
1871	!!gm bug ? : must be a do loop with mj0,mj1
1872
1873	! we duplicate factor scales for jj = 1 and jj = 2
1874	tg_e3t_0%d_value(:,1,:,1) = tg_e3t_0%d_value(:,2,:,1)
1875	tg_e3w_0%d_value(:,1,:,1) = tg_e3w_0%d_value(:,2,:,1)
1876	tg_e3u_0%d_value(:,1,:,1) = tg_e3u_0%d_value(:,2,:,1)
1877	tg_e3v_0%d_value(:,1,:,1) = tg_e3v_0%d_value(:,2,:,1)
1878	!tg_e3f_0%d_value(:,1,:,1) = tg_e3f_0%d_value(:,2,:,1)
1879
1880	! Control of the sign
1881	IF( MINVAL( tg_e3t_0%d_value (:,:,:,:) ) <= 0._dp ) CALL logger_fatal( ' GRID ZGR ZPS: e r r o r e3t_0 <= 0' )
1882	IF( MINVAL( tg_e3w_0%d_value (:,:,:,:) ) <= 0._dp ) CALL logger_fatal( ' GRID ZGR ZPS: e r r o r e3w_0 <= 0' )
1883	IF( MINVAL( tg_gdept_0%d_value(:,:,:,:) ) < 0._dp ) CALL logger_fatal( ' GRID ZGR ZPS: e r r o r gdept_0 < 0' )
1884	IF( MINVAL( tg_gdepw_0%d_value(:,:,:,:) ) < 0._dp ) CALL logger_fatal( ' GRID ZGR ZPS: e r r o r gdepw_0 < 0' )
1885
1886	!! Compute gdep3w_0 (vertical sum of e3w)
1887	!IF ( td_nam%l_isfcav ) THEN
1888	! ! if cavity
1889	! CALL grid_zgr__isf_fill_gdep3w_0(jpi, jpj, jpk, td_risfdep)
1890	!ELSE
1891	! ! no cavity
1892	! tg_gdep3w_0%d_value(:,:,1,1) = 0.5_dp * tg_e3w_0%d_value(:,:,1,1)
1893	! DO jk = 2, jpk
1894	! tg_gdep3w_0%d_value(:,:,jk,1) = tg_gdep3w_0%d_value(:,:,jk-1,1) + tg_e3w_0%d_value(:,:,jk,1)
1895	! END DO
1896	!END IF
1897
1898	END SUBROUTINE grid_zgr__zps_fill
1899	!-------------------------------------------------------------------
1900	!> @brief This subroutine check the bathymetry in levels
1901	!>
1902	!> @details
1903	!> ** Method : THe water column have to contained at least 2 cells
1904	!> Bathymetry and isfdraft are modified (dig/close) to respect
1905	!> this criterion.
1906	!>
1907	!>
1908	!> ** Action : - test compatibility between isfdraft and bathy
1909	!> - bathy and isfdraft are modified
1910	!>
1911	!> @author J.Paul
1912	!> @date September, 2015 - rewrite from zgr_isf
1913	!> @date October, 2016
1914	!> - add ice sheet coupling case
1915	!>
1916	!> @param[in] td_nam
1917	!> @param[in] jpi
1918	!> @param[in] jpj
1919	!> @param[in] jpk
1920	!> @param[in] td_bathy
1921	!> @param[in] td_risfdep
1922	!-------------------------------------------------------------------
1923	SUBROUTINE grid_zgr__isf_fill( td_nam, jpi,jpj,jpk, td_bathy, td_risfdep)
1924	IMPLICIT NONE
1925	! Argument
1926	TYPE(TNAMZ), INTENT(IN ) :: td_nam
1927	INTEGER(i4), INTENT(IN ) :: jpi
1928	INTEGER(i4), INTENT(IN ) :: jpj
1929	INTEGER(i4), INTENT(IN ) :: jpk
1930	TYPE(TVAR) , INTENT(INOUT) :: td_bathy
1931	TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
1932
1933	! local variable
1934	INTEGER(i4) :: ik, it
1935	INTEGER(i4) :: ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF)
1936
1937	INTEGER(i4), ALLOCATABLE, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)
1938
1939	REAL(dp) :: zdepth ! Ajusted ocean depth to avoid too small e3t
1940	REAL(dp) :: zmax ! Maximum and minimum depth
1941	REAL(dp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar
1942	REAL(dp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t
1943	REAL(dp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1944	REAL(dp) :: zdiff ! temporary scalar
1945
1946	REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: zrisfdep, zmask ! 2D workspace (ISH)
1947
1948	! loop indices
1949	INTEGER(i4) :: ji
1950	INTEGER(i4) :: jj
1951	INTEGER(i4) :: jk
1952	INTEGER(i4) :: jl
1953
1954	!----------------------------------------------------------------
1955
1956	! (ISF) compute misfdep
1957	WHERE( td_risfdep%d_value(:,:,1,1) == 0._dp .AND. &
1958	& td_bathy%d_value (:,:,1,1) /= 0 )
1959	! open water : set misfdep to 1
1960	tg_misfdep%d_value(:,:,1,1) = 1
1961	ELSEWHERE
1962	! iceshelf : initialize misfdep to second level
1963	tg_misfdep%d_value(:,:,1,1) = 2
1964	END WHERE
1965
1966	ALLOCATE(zmask (jpi,jpj))
1967	ALLOCATE(zrisfdep(jpi,jpj))
1968	ALLOCATE(zmisfdep(jpi,jpj))
1969
1970	! Compute misfdep for ocean points (i.e. first wet level)
1971	! find the first ocean level such that the first level thickness
1972	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
1973	! e3t_0 is the reference level thickness
1974	DO jk = 2, jpk-1
1975	zdepth = tg_gdepw_1d%d_value(1,1,jk+1,1) - MIN( td_nam%d_e3zps_min, &
1976	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,jk,1) )
1977	WHERE( 0._dp < td_risfdep%d_value(:,:,1,1) .AND. &
1978	& td_risfdep%d_value(:,:,1,1) >= zdepth )
1979	tg_misfdep%d_value(:,:,1,1) = jk+1
1980	END WHERE
1981	END DO
1982
1983	WHERE( 0._dp < td_risfdep%d_value(:,:,1,1) .AND. &
1984	& td_risfdep%d_value(:,:,1,1) <= tg_e3t_1d%d_value(1,1,1,1) )
1985	td_risfdep%d_value(:,:,1,1) = 0.
1986	tg_misfdep%d_value(:,:,1,1) = 1
1987	END WHERE
1988
1989	! remove very shallow ice shelf (less than ~ 10m if 75L)
1990	WHERE( td_risfdep%d_value(:,:,1,1) <= 10._dp .AND. &
1991	& tg_misfdep%d_value(:,:,1,1) > 1 )
1992	tg_misfdep%d_value(:,:,1,1) = 0
1993	td_risfdep%d_value(:,:,1,1) = 0.0_dp
1994	tg_mbathy%d_value (:,:,1,1) = 0
1995	td_bathy%d_value (:,:,1,1) = 0.0_dp
1996	END WHERE
1997	WHERE( td_bathy%d_value(:,:,1,1) <= 30.0_dp .AND. &
1998	& tg_gphit%d_value(:,:,1,1) < -60._dp )
1999	tg_misfdep%d_value(:,:,1,1) = 0
2000	td_risfdep%d_value(:,:,1,1) = 0.0_dp
2001	tg_mbathy%d_value (:,:,1,1) = 0
2002	td_bathy%d_value (:,:,1,1) = 0.0_dp
2003	END WHERE
2004
2005	! basic check for the compatibility of bathy and risfdep.
2006	! I think it should be offline because it is not perfect and cannot solved all the situation
2007	! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
2008	DO jl = 1, 10
2009
2010	WHERE( td_bathy%d_value(:,:,1,1) <= td_risfdep%d_value(:,:,1,1) + td_nam%d_isfhmin )
2011	tg_misfdep%d_value(:,:,1,1) = 0
2012	td_risfdep%d_value(:,:,1,1) = 0._dp
2013	tg_mbathy%d_value (:,:,1,1) = 0
2014	td_bathy%d_value (:,:,1,1) = 0._dp
2015	END WHERE
2016
2017	WHERE( tg_mbathy%d_value(:,:,1,1) <= 0 )
2018	tg_misfdep%d_value(:,:,1,1) = 0
2019	td_risfdep%d_value(:,:,1,1) = 0._dp
2020	tg_mbathy%d_value (:,:,1,1) = 0
2021	td_bathy%d_value (:,:,1,1) = 0._dp
2022	ENDWHERE
2023
2024	!! lk_mpp not added
2025
2026	IF( td_nam%i_perio == 1 .OR. &
2027	& td_nam%i_perio == 4 .OR. &
2028	& td_nam%i_perio == 6 )THEN
2029	! local domain is cyclic east-west
2030	tg_misfdep%d_value( 1 ,:,1,1) = tg_misfdep%d_value(jpi-1,:,1,1)
2031	tg_misfdep%d_value(jpi,:,1,1) = tg_misfdep%d_value( 2 ,:,1,1)
2032
2033	tg_mbathy%d_value( 1 ,:,1,1) = tg_mbathy%d_value(jpi-1,:,1,1)
2034	tg_mbathy%d_value(jpi,:,1,1) = tg_mbathy%d_value( 2 ,:,1,1)
2035	ENDIF
2036
2037	! split last cell if possible (only where water column is 2 cell or less)
2038	! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
2039	IF( .NOT. td_nam%l_iscpl )THEN
2040	DO jk = jpk-1, 1, -1
2041	zmax = tg_gdepw_1d%d_value(1,1,jk,1) + &
2042	& MIN( td_nam%d_e3zps_min, &
2043	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,jk,1) )
2044
2045	WHERE( td_bathy%d_value(:,:,1,1) > tg_gdepw_1d%d_value(1,1,jk,1).AND. &
2046	& td_bathy%d_value(:,:,1,1) <= zmax .AND. &
2047	& tg_misfdep%d_value(:,:,1,1) + 1 >= tg_mbathy%d_value(:,:,1,1) )
2048
2049	tg_mbathy%d_value(:,:,1,1) = jk
2050	td_bathy%d_value (:,:,1,1) = zmax
2051
2052	END WHERE
2053	END DO
2054	ENDIF
2055
2056	! split top cell if possible (only where water column is 2 cell or less)
2057	DO jk = 2, jpk-1
2058	zmax = tg_gdepw_1d%d_value(1,1,jk+1,1) - &
2059	& MIN( td_nam%d_e3zps_min, &
2060	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,jk,1) )
2061
2062	WHERE( td_risfdep%d_value(:,:,1,1) < tg_gdepw_1d%d_value(1,1,jk+1,1) .AND. &
2063	& td_risfdep%d_value(:,:,1,1) >= zmax .AND. &
2064	& tg_misfdep%d_value(:,:,1,1) + 1 >= tg_mbathy%d_value(:,:,1,1) )
2065
2066	tg_misfdep%d_value(:,:,1,1) = jk
2067	td_risfdep%d_value(:,:,1,1) = zmax
2068
2069	END WHERE
2070	END DO
2071
2072	! Case where bathy and risfdep compatible but not the level
2073	! variable mbathy/misfdep because of partial cell condition
2074	DO jj = 1, jpj
2075	DO ji = 1, jpi
2076	! find the minimum change option:
2077	! test bathy
2078	IF( td_risfdep%d_value(ji,jj,1,1) > 1 )THEN
2079	IF( .NOT. td_nam%l_iscpl )THEN
2080
2081	ik=tg_mbathy%d_value(ji,jj,1,1)
2082	zbathydiff = ABS( td_bathy%d_value(ji,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2083	& + MIN(td_nam%d_e3zps_min, &
2084	& td_nam%d_e3zps_rat * tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2085
2086	ik=tg_misfdep%d_value(ji,jj,1,1)
2087	zrisfdepdiff = ABS( td_risfdep%d_value(ji,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2088	& - MIN( td_nam%d_e3zps_min, &
2089	& td_nam%d_e3zps_rat * tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2090
2091	IF( td_bathy%d_value (ji,jj,1,1) > td_risfdep%d_value (ji,jj,1,1) .AND. &
2092	& tg_mbathy%d_value(ji,jj,1,1) < tg_misfdep%d_value(ji,jj,1,1) )THEN
2093
2094	IF( zbathydiff <= zrisfdepdiff )THEN
2095
2096	tg_mbathy%d_value(ji,jj,1,1) = tg_mbathy%d_value(ji,jj,1,1) + 1
2097
2098	ik=tg_mbathy%d_value(ji,jj,1,1)
2099	td_bathy%d_value (ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2100	& MIN( td_nam%d_e3zps_min, &
2101	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2102	ELSE
2103
2104	tg_misfdep%d_value(ji,jj,1,1) = tg_misfdep%d_value(ji,jj,1,1) - 1
2105
2106	ik=tg_misfdep%d_value(ji,jj,1,1)
2107	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) - &
2108	& MIN( td_nam%d_e3zps_min, &
2109	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1) )
2110	ENDIF
2111
2112	ENDIF
2113
2114	ELSE
2115
2116	IF( td_bathy%d_value (ji,jj,1,1) > td_risfdep%d_value (ji,jj,1,1) .AND. &
2117	& tg_mbathy%d_value(ji,jj,1,1) < tg_misfdep%d_value(ji,jj,1,1) )THEN
2118
2119	tg_misfdep%d_value(ji,jj,1,1) = tg_misfdep%d_value(ji,jj,1,1) - 1
2120
2121	ik=tg_misfdep%d_value(ji,jj,1,1)
2122	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) - &
2123	& MIN( td_nam%d_e3zps_min, &
2124	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1) )
2125	ENDIF
2126
2127	ENDIF
2128
2129	ENDIF
2130	ENDDO
2131	ENDDO
2132
2133	! At least 2 levels for water thickness at T, U, and V point.
2134	DO jj = 1, jpj
2135	DO ji = 1, jpi
2136	! find the minimum change option:
2137	! test bathy
2138	IF( tg_misfdep%d_value(ji,jj,1,1) == tg_mbathy%d_value(ji,jj,1,1) .AND. &
2139	& tg_mbathy%d_value (ji,jj,1,1) > 1) THEN
2140
2141	IF ( .NOT. td_nam%l_iscpl ) THEN
2142
2143	ik=tg_mbathy%d_value(ji,jj,1,1)
2144	zbathydiff = ABS( td_bathy%d_value(ji,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2145	& + MIN( td_nam%d_e3zps_min, &
2146	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2147
2148	ik=tg_misfdep%d_value(ji,jj,1,1)
2149	zrisfdepdiff = ABS( td_risfdep%d_value(ji,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2150	& - MIN( td_nam%d_e3zps_min, &
2151	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2152
2153	IF( zbathydiff <= zrisfdepdiff )THEN
2154
2155	tg_mbathy%d_value(ji,jj,1,1) = tg_mbathy%d_value(ji,jj,1,1) + 1
2156
2157	ik=tg_mbathy%d_value(ji,jj,1,1)
2158	td_bathy%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2159	& MIN( td_nam%d_e3zps_min, &
2160	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2161	ELSE
2162
2163	tg_misfdep%d_value(ji,jj,1,1) = tg_misfdep%d_value(ji,jj,1,1) - 1
2164
2165	ik=tg_misfdep%d_value(ji,jj,1,1)
2166	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2167	& MIN( td_nam%d_e3zps_min, &
2168	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2169	ENDIF
2170
2171	ELSE
2172
2173	tg_misfdep%d_value(ji,jj,1,1) = tg_misfdep%d_value(ji,jj,1,1) - 1
2174
2175	ik=tg_misfdep%d_value(ji,jj,1,1)
2176	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2177	& MIN( td_nam%d_e3zps_min, &
2178	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2179
2180	ENDIF
2181	ENDIF
2182
2183	ENDDO
2184	ENDDO
2185
2186	! point V mbathy(ji,jj ) == misfdep(ji,jj+1)
2187	DO jj = 1, jpj-1
2188	DO ji = 1, jpi-1
2189
2190	IF( tg_mbathy%d_value(ji,jj ,1,1) == tg_misfdep%d_value(ji,jj+1,1,1) .AND. &
2191	& tg_mbathy%d_value(ji,jj ,1,1) > 1 )THEN
2192	IF ( .NOT. td_nam%l_iscpl ) THEN
2193
2194	ik=tg_mbathy%d_value(ji,jj ,1,1)
2195	zbathydiff = ABS( td_bathy%d_value(ji,jj ,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2196	& + MIN( td_nam%d_e3zps_min, &
2197	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2198
2199	ik=tg_misfdep%d_value(ji,jj+1,1,1)
2200	zrisfdepdiff = ABS( td_risfdep%d_value(ji,jj+1,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2201	& - MIN( td_nam%d_e3zps_min, &
2202	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2203
2204	IF( zbathydiff <= zrisfdepdiff )THEN
2205
2206	tg_mbathy%d_value(ji,jj ,1,1) = tg_mbathy%d_value(ji,jj ,1,1) + 1
2207
2208	ik=tg_mbathy%d_value(ji,jj ,1,1)
2209	td_bathy%d_value (ji,jj ,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2210	& MIN( td_nam%d_e3zps_min, &
2211	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2212	ELSE
2213
2214	tg_misfdep%d_value(ji,jj+1,1,1) = tg_misfdep%d_value(ji,jj+1,1,1) - 1
2215
2216	ik=tg_misfdep%d_value(ji,jj+1,1,1)
2217	td_risfdep%d_value(ji,jj+1,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2218	& MIN( td_nam%d_e3zps_min, &
2219	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2220	ENDIF
2221
2222	ELSE
2223	tg_misfdep%d_value(ji,jj+1,1,1) = tg_misfdep%d_value(ji,jj+1,1,1) - 1
2224
2225	ik=tg_misfdep%d_value(ji,jj+1,1,1)
2226	td_risfdep%d_value(ji,jj+1,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2227	& MIN( td_nam%d_e3zps_min, &
2228	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2229	ENDIF
2230
2231	ENDIF
2232
2233	ENDDO
2234	ENDDO
2235
2236	!! lk_mpp not added
2237
2238	! point V mbathy(ji,jj+1) == misfdep(ji,jj )
2239	DO jj = 1, jpj-1
2240	DO ji = 1, jpi-1
2241
2242	IF( tg_mbathy%d_value(ji,jj+1,1,1) == tg_misfdep%d_value(ji,jj ,1,1) .AND. &
2243	& tg_mbathy%d_value(ji,jj ,1,1) > 1) THEN
2244	IF ( .NOT. td_nam%l_iscpl ) THEN
2245
2246	ik=tg_mbathy%d_value(ji,jj+1,1,1)
2247	zbathydiff = ABS( td_bathy%d_value(ji,jj+1,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2248	& + MIN( td_nam%d_e3zps_min, &
2249	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2250
2251	ik=tg_misfdep%d_value(ji,jj ,1,1)
2252	zrisfdepdiff = ABS( td_risfdep%d_value(ji,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2253	& - MIN( td_nam%d_e3zps_min, &
2254	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2255
2256	IF( zbathydiff <= zrisfdepdiff )THEN
2257
2258	tg_mbathy%d_value (ji,jj+1,1,1) = tg_mbathy%d_value(ji,jj+1,1,1) + 1
2259
2260	ik=tg_mbathy%d_value(ji,jj+1,1,1)
2261	td_bathy%d_value (ji,jj+1,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2262	& MIN( td_nam%d_e3zps_min, &
2263	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2264	ELSE
2265
2266	tg_misfdep%d_value(ji,jj ,1,1) = tg_misfdep%d_value(ji,jj ,1,1) - 1
2267
2268	ik=tg_misfdep%d_value(ji,jj ,1,1)
2269	td_risfdep%d_value(ji,jj ,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2270	& MIN( td_nam%d_e3zps_min, &
2271	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2272	END IF
2273
2274	ELSE
2275
2276	tg_misfdep%d_value(ji,jj ,1,1) = tg_misfdep%d_value(ji,jj ,1,1) - 1
2277
2278	ik=tg_misfdep%d_value(ji,jj ,1,1)
2279	td_risfdep%d_value(ji,jj ,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2280	& MIN( td_nam%d_e3zps_min, &
2281	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2282
2283	ENDIF
2284
2285	ENDIF
2286
2287	ENDDO
2288	ENDDO
2289
2290	!! lk_mpp not added
2291
2292	! point U mbathy(ji ,jj) EQ misfdep(ji+1,jj)
2293	DO jj = 1, jpj-1
2294	DO ji = 1, jpi-1
2295
2296	IF( tg_mbathy%d_value(ji ,jj,1,1) == tg_misfdep%d_value(ji+1,jj,1,1) .AND. &
2297	& tg_mbathy%d_value(ji ,jj,1,1) > 1 )THEN
2298	IF ( .NOT. td_nam%l_iscpl ) THEN
2299
2300	ik=tg_mbathy%d_value(ji ,jj,1,1)
2301	zbathydiff = ABS( td_bathy%d_value(ji ,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2302	& + MIN( td_nam%d_e3zps_min, &
2303	& td_nam%d_e3zps_rat* tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2304
2305	ik=tg_misfdep%d_value(ji+1,jj,1,1)
2306	zrisfdepdiff = ABS( td_risfdep%d_value(ji+1,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2307	& - MIN( td_nam%d_e3zps_min, &
2308	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2309
2310	IF( zbathydiff <= zrisfdepdiff )THEN
2311
2312	tg_mbathy%d_value (ji ,jj,1,1) = tg_mbathy%d_value(ji ,jj,1,1) + 1
2313
2314	ik=tg_mbathy%d_value(ji ,jj,1,1)
2315	td_bathy%d_value (ji ,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2316	& MIN( td_nam%d_e3zps_min, &
2317	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2318	ELSE
2319	tg_misfdep%d_value(ji+1,jj,1,1) = tg_misfdep%d_value(ji+1,jj,1,1) - 1
2320
2321	ik=tg_misfdep%d_value(ji+1,jj,1,1)
2322	td_risfdep%d_value(ji+1,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2323	& MIN( td_nam%d_e3zps_min, &
2324	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2325	END IF
2326
2327	ELSE
2328	tg_misfdep%d_value(ji+1,jj,1,1)= tg_misfdep%d_value(ji+1,jj,1,1) - 1
2329
2330	ik=tg_misfdep%d_value(ji+1,jj,1,1)
2331	td_risfdep%d_value(ji+1,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2332	& MIN( td_nam%d_e3zps_min, &
2333	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2334	ENDIF
2335	ENDIF
2336
2337	ENDDO
2338	ENDDO
2339
2340	!! lk_mpp not added
2341
2342	! point U mbathy(ji+1,jj) == misfdep(ji ,jj)
2343	DO jj = 1, jpj-1
2344	DO ji = 1, jpi-1
2345
2346	IF( tg_mbathy%d_value(ji+1,jj,1,1) == tg_misfdep%d_value(ji ,jj,1,1) .AND. &
2347	& tg_mbathy%d_value(ji ,jj,1,1) > 1 )THEN
2348	IF ( .NOT. td_nam%l_iscpl ) THEN
2349
2350	ik=tg_mbathy%d_value(ji+1,jj,1,1)
2351	zbathydiff = ABS( td_bathy%d_value(ji+1,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik+1,1) &
2352	& + MIN( td_nam%d_e3zps_min, &
2353	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1)) ))
2354
2355	ik=tg_misfdep%d_value(ji ,jj,1,1)
2356	zrisfdepdiff = ABS( td_risfdep%d_value(ji ,jj,1,1) - ( tg_gdepw_1d%d_value(1,1,ik,1) &
2357	& - MIN( td_nam%d_e3zps_min, &
2358	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik-1,1)) ))
2359
2360	IF( zbathydiff <= zrisfdepdiff )THEN
2361
2362	tg_mbathy%d_value (ji+1,jj,1,1) = tg_mbathy%d_value(ji+1,jj,1,1) + 1
2363
2364	ik=tg_mbathy%d_value(ji+1,jj,1,1)
2365	td_bathy%d_value (ji+1,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2366	& MIN( td_nam%d_e3zps_min, &
2367	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik+1,1) )
2368	ELSE
2369	tg_misfdep%d_value(ji ,jj,1,1) = tg_misfdep%d_value(ji ,jj,1,1) - 1
2370
2371	ik=tg_misfdep%d_value(ji ,jj,1,1)
2372	td_risfdep%d_value(ji ,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2373	& MIN( td_nam%d_e3zps_min, &
2374	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2375	END IF
2376
2377	ELSE
2378
2379	tg_misfdep%d_value(ji ,jj,1,1) = tg_misfdep%d_value(ji ,jj,1,1) - 1
2380
2381	ik=tg_misfdep%d_value(ji ,jj,1,1)
2382	td_risfdep%d_value(ji ,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2383	& MIN( td_nam%d_e3zps_min, &
2384	& td_nam%d_e3zps_rat*tg_e3t_1d%d_value(1,1,ik,1) )
2385
2386	ENDIF
2387
2388	ENDIF
2389
2390	ENDDO
2391	ENDDO
2392
2393	END DO ! jl
2394	! end dig bathy/ice shelf to be compatible
2395
2396	! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
2397	DO jl = 1,20
2398
2399	! remove single point "bay" on isf coast line in the ice shelf draft'
2400	DO jk = 2, jpk
2401	WHERE( tg_misfdep%d_value(:,:,1,1) == 0 )
2402	tg_misfdep%d_value(:,:,1,1)=jpk
2403	END WHERE
2404
2405	zmask(:,:)=0
2406	WHERE( tg_misfdep%d_value(:,:,1,1) <= jk ) zmask(:,:)=1
2407
2408	DO jj = 2, jpj-1
2409	DO ji = 2, jpi-1
2410	IF( tg_misfdep%d_value(ji,jj,1,1) == jk )THEN
2411
2412	ibtest = zmask(ji-1,jj ) &
2413	& + zmask(ji+1,jj ) &
2414	& + zmask(ji ,jj-1) &
2415	& + zmask(ji ,jj+1)
2416	IF( ibtest <= 1 )THEN
2417	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,jk+1,1)
2418	tg_misfdep%d_value(ji,jj,1,1) = jk+1
2419	IF( tg_misfdep%d_value(ji,jj,1,1) > tg_mbathy%d_value(ji,jj,1,1) )THEN
2420	tg_misfdep%d_value(ji,jj,1,1) = jpk
2421	ENDIF
2422	ENDIF
2423
2424	ENDIF
2425	ENDDO
2426	ENDDO
2427	ENDDO
2428
2429	WHERE( tg_misfdep%d_value(:,:,1,1) == jpk )
2430	tg_misfdep%d_value(:,:,1,1)=0
2431	td_risfdep%d_value(:,:,1,1)=0.
2432	tg_mbathy%d_value (:,:,1,1)=0
2433	td_bathy%d_value (:,:,1,1)=0.
2434	END WHERE
2435
2436	!! lk_mpp not added
2437
2438	! remove single point "bay" on bathy coast line beneath an ice shelf'
2439	DO jk = jpk,1,-1
2440
2441	zmask(:,:)=0
2442	WHERE( tg_mbathy%d_value(:,:,1,1) >= jk ) zmask(:,:)=1
2443
2444	DO jj = 2, jpj-1
2445	DO ji = 2, jpi-1
2446	IF( tg_mbathy%d_value (ji,jj,1,1) == jk .AND. &
2447	& tg_misfdep%d_value(ji,jj,1,1) >= 2 )THEN
2448
2449	ibtest = zmask(ji-1,jj ) &
2450	& + zmask(ji+1,jj ) &
2451	& + zmask(ji ,jj-1) &
2452	& + zmask(ji ,jj+1)
2453	IF( ibtest <= 1 )THEN
2454	td_bathy%d_value (ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,jk,1)
2455	tg_mbathy%d_value(ji,jj,1,1) = jk-1
2456	IF( tg_misfdep%d_value(ji,jj,1,1) > tg_mbathy%d_value(ji,jj,1,1) )THEN
2457	tg_mbathy%d_value(ji,jj,1,1) = 0
2458	ENDIF
2459	ENDIF
2460
2461	ENDIF
2462	ENDDO
2463	ENDDO
2464	ENDDO
2465
2466	WHERE( tg_mbathy%d_value(:,:,1,1) == 0 )
2467	tg_misfdep%d_value(:,:,1,1)=0
2468	td_risfdep%d_value(:,:,1,1)=0.
2469	tg_mbathy%d_value (:,:,1,1)=0
2470	td_bathy%d_value (:,:,1,1)=0.
2471	END WHERE
2472
2473	!! lk_mpp not added
2474
2475	! fill hole in ice shelf
2476	zmisfdep(:,:) = tg_misfdep%d_value(:,:,1,1)
2477	zrisfdep(:,:) = td_risfdep%d_value(:,:,1,1)
2478	WHERE( zmisfdep(:,:) <= 1 ) zmisfdep(:,:)=jpk
2479
2480	DO jj = 2, jpj-1
2481	DO ji = 2, jpi-1
2482
2483	ibtestim1 = zmisfdep(ji-1,jj )
2484	ibtestip1 = zmisfdep(ji+1,jj )
2485
2486	ibtestjm1 = zmisfdep(ji ,jj-1)
2487	ibtestjp1 = zmisfdep(ji ,jj+1)
2488
2489	IF( zmisfdep(ji,jj) >= tg_mbathy%d_value(ji-1,jj ,1,1) ) ibtestim1 = jpk
2490	IF( zmisfdep(ji,jj) >= tg_mbathy%d_value(ji+1,jj ,1,1) ) ibtestip1 = jpk
2491	IF( zmisfdep(ji,jj) >= tg_mbathy%d_value(ji ,jj-1,1,1) ) ibtestjm1 = jpk
2492	IF( zmisfdep(ji,jj) >= tg_mbathy%d_value(ji ,jj+1,1,1) ) ibtestjp1 = jpk
2493
2494	ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
2495	IF( ibtest == jpk .AND. &
2496	& tg_misfdep%d_value(ji,jj,1,1) >= 2 )THEN
2497	tg_mbathy%d_value (ji,jj,1,1) = 0
2498	td_bathy%d_value (ji,jj,1,1) = 0.0_dp
2499	tg_misfdep%d_value(ji,jj,1,1) = 0
2500	td_risfdep%d_value(ji,jj,1,1) = 0.0_dp
2501	END IF
2502
2503	IF( zmisfdep(ji,jj) < ibtest .AND. &
2504	& tg_misfdep%d_value(ji,jj,1,1) >= 2 )THEN
2505	tg_misfdep%d_value(ji,jj,1,1) = ibtest
2506	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ibtest,1)
2507	ENDIF
2508
2509	ENDDO
2510	ENDDO
2511
2512	!! lk_mpp not added
2513
2514	!! fill hole in bathymetry
2515	zmbathy(:,:) = tg_mbathy%d_value(:,:,1,1)
2516	DO jj = 2, jpj-1
2517	DO ji = 2, jpi-1
2518
2519	ibtestim1 = zmbathy(ji-1,jj )
2520	ibtestip1 = zmbathy(ji+1,jj )
2521
2522	ibtestjm1 = zmbathy(ji ,jj-1)
2523	ibtestjp1 = zmbathy(ji ,jj+1)
2524
2525	IF( zmbathy(ji,jj) < tg_misfdep%d_value(ji-1,jj ,1,1) ) ibtestim1 = 0
2526	IF( zmbathy(ji,jj) < tg_misfdep%d_value(ji+1,jj ,1,1) ) ibtestip1 = 0
2527	IF( zmbathy(ji,jj) < tg_misfdep%d_value(ji ,jj-1,1,1) ) ibtestjm1 = 0
2528	IF( zmbathy(ji,jj) < tg_misfdep%d_value(ji ,jj+1,1,1) ) ibtestjp1 = 0
2529
2530	ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
2531	IF( ibtest == 0 .AND. tg_misfdep%d_value(ji,jj,1,1) >= 2) THEN
2532	tg_mbathy%d_value (ji,jj,1,1) = 0
2533	td_bathy%d_value (ji,jj,1,1) = 0.0_dp
2534	tg_misfdep%d_value(ji,jj,1,1) = 0
2535	td_risfdep%d_value(ji,jj,1,1) = 0.0_dp
2536	END IF
2537
2538	IF( ibtest < zmbathy(ji,jj) .AND. &
2539	& tg_misfdep%d_value(ji,jj,1,1) >= 2) THEN
2540	tg_mbathy%d_value(ji,jj,1,1) = ibtest
2541	td_bathy%d_value (ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ibtest+1,1)
2542	ENDIF
2543
2544	ENDDO
2545	ENDDO
2546
2547	!! lk_mpp not added
2548
2549	! if not compatible after all check (ie U point water column less than 2 cells), mask U
2550	DO jj = 1, jpj-1
2551	DO ji = 1, jpi-1
2552	IF( tg_mbathy%d_value(ji ,jj,1,1) == tg_misfdep%d_value(ji+1,jj,1,1) .AND. &
2553	& tg_mbathy%d_value(ji ,jj,1,1) >= 1 .AND. &
2554	& tg_mbathy%d_value(ji+1,jj,1,1) >= 1 )THEN
2555
2556	tg_mbathy%d_value(ji,jj,1,1) = tg_mbathy%d_value(ji,jj,1,1) - 1
2557
2558	ik=tg_mbathy%d_value(ji,jj,1,1)
2559	td_bathy%d_value (ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
2560	ENDIF
2561	ENDDO
2562	ENDDO
2563
2564	!! lk_mpp not added
2565
2566	! if not compatible after all check (ie U point water column less than 2 cells), mask U
2567	DO jj = 1, jpj-1
2568	DO ji = 1, jpi-1
2569	IF( tg_misfdep%d_value(ji ,jj,1,1) == tg_mbathy%d_value(ji+1,jj,1,1) .AND. &
2570	& tg_mbathy%d_value (ji ,jj,1,1) >= 1 .AND.&
2571	& tg_mbathy%d_value (ji+1,jj,1,1) >= 1 )THEN
2572
2573	tg_mbathy%d_value(ji+1,jj,1,1) = tg_mbathy%d_value(ji+1,jj,1,1) - 1
2574
2575	ik=tg_mbathy%d_value(ji+1,jj,1,1)
2576	td_bathy%d_value(ji+1,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
2577	ENDIF
2578	ENDDO
2579	ENDDO
2580
2581	!! lk_mpp not added
2582
2583	! if not compatible after all check (ie V point water column less than 2 cells), mask V
2584	DO jj = 1, jpj-1
2585	DO ji = 1, jpi
2586	IF( tg_mbathy%d_value(ji,jj ,1,1) == tg_misfdep%d_value(ji,jj+1,1,1) .AND. &
2587	& tg_mbathy%d_value(ji,jj ,1,1) >= 1 .AND. &
2588	& tg_mbathy%d_value(ji,jj+1,1,1) >= 1 )THEN
2589
2590	tg_mbathy%d_value(ji,jj,1,1) = tg_mbathy%d_value(ji,jj,1,1) - 1
2591
2592	ik=tg_mbathy%d_value(ji,jj,1,1)
2593	td_bathy%d_value (ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
2594
2595	ENDIF
2596	ENDDO
2597	ENDDO
2598
2599	!! lk_mpp not added
2600
2601	! if not compatible after all check (ie V point water column less than 2 cells), mask V
2602	DO jj = 1, jpj-1
2603	DO ji = 1, jpi
2604	IF( tg_misfdep%d_value(ji,jj ,1,1) == tg_mbathy%d_value(ji,jj+1,1,1) .AND.&
2605	& tg_mbathy%d_value (ji,jj ,1,1) >= 1 .AND.&
2606	& tg_mbathy%d_value (ji,jj+1,1,1) >= 1 )THEN
2607
2608	tg_mbathy%d_value(ji,jj+1,1,1) = tg_mbathy%d_value(ji,jj+1,1,1) - 1
2609
2610	ik=tg_mbathy%d_value(ji,jj+1,1,1)
2611	td_bathy%d_value (ji,jj+1,1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
2612	ENDIF
2613	ENDDO
2614	ENDDO
2615
2616	!! lk_mpp not added
2617
2618	! if not compatible after all check, mask T
2619	DO jj = 1, jpj
2620	DO ji = 1, jpi
2621	IF( tg_mbathy%d_value(ji,jj,1,1) <= tg_misfdep%d_value(ji,jj,1,1) )THEN
2622	tg_misfdep%d_value(ji,jj,1,1) = 0
2623	td_risfdep%d_value(ji,jj,1,1) = 0._dp
2624	tg_mbathy%d_value (ji,jj,1,1) = 0
2625	td_bathy%d_value (ji,jj,1,1) = 0._dp
2626	ENDIF
2627	ENDDO
2628	ENDDO
2629
2630	WHERE( tg_mbathy%d_value(:,:,1,1) == 1 )
2631	tg_mbathy%d_value (:,:,1,1) = 0
2632	td_bathy%d_value (:,:,1,1) = 0.0_dp
2633	tg_misfdep%d_value(:,:,1,1) = 0
2634	td_risfdep%d_value(:,:,1,1) = 0.0_dp
2635	END WHERE
2636	ENDDO
2637	! end check compatibility ice shelf/bathy
2638
2639	! remove very shallow ice shelf (less than ~ 10m if 75L)
2640	WHERE( td_risfdep%d_value(:,:,1,1) <= 10._dp )
2641	tg_misfdep%d_value(:,:,1,1) = 1
2642	td_risfdep%d_value(:,:,1,1) = 0.0_dp;
2643	END WHERE
2644
2645	DEALLOCATE(zmask )
2646	DEALLOCATE(zrisfdep)
2647	DEALLOCATE(zmisfdep)
2648
2649	! compute scale factor and depth at T- and W- points
2650	DO jj = 1, jpj
2651	DO ji = 1, jpi
2652	ik = tg_mbathy%d_value(ji,jj,1,1)
2653	IF( ik > 0 ) THEN ! ocean point only
2654	! max ocean level case
2655	IF( ik == jpk-1 ) THEN
2656
2657	zdepwp = td_bathy%d_value(ji,jj,1,1)
2658	ze3tp = td_bathy%d_value(ji,jj,1,1) - tg_gdepw_1d%d_value(1,1,ik,1)
2659	ze3wp = 0.5_dp * tg_e3w_1d%d_value(1,1,ik,1) &
2660	& * ( 1._dp + ( ze3tp/tg_e3t_1d%d_value(1,1,ik,1) ) )
2661	tg_e3t_0%d_value (ji,jj,ik ,1) = ze3tp
2662	tg_e3t_0%d_value (ji,jj,ik+1,1) = ze3tp
2663	tg_e3w_0%d_value (ji,jj,ik ,1) = ze3wp
2664	tg_e3w_0%d_value (ji,jj,ik+1,1) = ze3wp
2665
2666	tg_gdepw_0%d_value(ji,jj,ik+1,1) = zdepwp
2667	tg_gdept_0%d_value(ji,jj,ik ,1) = tg_gdept_1d%d_value(1 ,1 ,ik-1,1) + ze3wp
2668	tg_gdept_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value (ji,jj,ik ,1) + ze3tp
2669	!
2670	ELSE ! standard case
2671
2672	IF( td_bathy%d_value(ji,jj,1,1) <= tg_gdepw_1d%d_value(1,1,ik+1,1) ) THEN
2673	tg_gdepw_0%d_value(ji,jj,ik+1,1) = td_bathy%d_value(ji,jj,1,1)
2674	ELSE
2675	tg_gdepw_0%d_value(ji,jj,ik+1,1) = tg_gdepw_1d%d_value(1,1,ik+1,1)
2676	ENDIF
2677	!
2678	!gm Bug? check the gdepw_1d
2679	! ... on ik
2680	tg_gdept_0%d_value(ji,jj,ik,1) = tg_gdepw_1d%d_value(1,1,ik,1) + &
2681	& ( tg_gdepw_0%d_value(ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
2682	& * ( ( tg_gdept_1d%d_value(1,1, ik ,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) &
2683	& / ( tg_gdepw_1d%d_value(1,1, ik+1,1) - tg_gdepw_1d%d_value(1,1,ik,1) ) )
2684
2685	tg_e3t_0%d_value(ji,jj,ik ,1) = tg_gdepw_0%d_value(ji,jj,ik+1,1) - tg_gdepw_1d%d_value(1,1,ik ,1)
2686	tg_e3w_0%d_value(ji,jj,ik ,1) = tg_gdept_0%d_value(ji,jj,ik ,1) - tg_gdept_1d%d_value(1,1,ik-1,1)
2687	! ... on ik+1
2688	tg_e3w_0%d_value(ji,jj,ik+1,1) = tg_e3t_0%d_value(ji,jj,ik,1)
2689	tg_e3t_0%d_value(ji,jj,ik+1,1) = tg_e3t_0%d_value(ji,jj,ik,1)
2690	!
2691	ENDIF
2692	ENDIF
2693	END DO
2694	END DO
2695	!
2696	it = 0
2697	DO jj = 1, jpj
2698	DO ji = 1, jpi
2699	ik = tg_mbathy%d_value(ji,jj,1,1)
2700	IF( ik > 0 ) THEN ! ocean point only
2701	tg_e3tp%d_value(ji,jj,1,1) = tg_e3t_0%d_value(ji,jj,ik,1)
2702	tg_e3wp%d_value(ji,jj,1,1) = tg_e3w_0%d_value(ji,jj,ik,1)
2703	! test
2704	zdiff= tg_gdepw_0%d_value(ji,jj,ik+1,1) - tg_gdept_0%d_value(ji,jj,ik,1)
2705	IF( zdiff <= 0._dp ) THEN
2706	it = it + 1
2707	CALL logger_info(' it = '//TRIM(fct_str(it))//&
2708	& ' ik = '//TRIM(fct_str(ik))//&
2709	& ' (i,j) = '//trim(fct_str(ji))//' '//TRIM(fct_str(jj)))
2710	CALL logger_info(' bathy = '//TRIM(fct_str(td_bathy%d_value(ji,jj,1,1))))
2711	CALL logger_info(' gdept_0 = '//TRIM(fct_str(tg_gdept_0%d_value(ji,jj,ik,1)))//&
2712	& ' gdepw_0 = '//TRIM(fct_str(tg_gdepw_0%d_value(ji,jj,ik+1,1)))//&
2713	& ' zdiff = '//TRIM(fct_str(zdiff)))
2714	CALL logger_info(' e3tp = '//TRIM(fct_str(tg_e3t_0%d_value(ji,jj,ik,1)))//&
2715	& ' e3wp = '//TRIM(fct_str(tg_e3w_0%d_value(ji,jj,ik,1))))
2716	ENDIF
2717	ENDIF
2718	END DO
2719	END DO
2720	!
2721	! (ISF) Definition of e3t, u, v, w for ISF case
2722	DO jj = 1, jpj
2723	DO ji = 1, jpi
2724	ik = tg_misfdep%d_value(ji,jj,1,1)
2725	IF( ik > 1 ) THEN ! ice shelf point only
2726
2727	IF( td_risfdep%d_value(ji,jj,1,1) < tg_gdepw_1d%d_value(1,1,ik,1) )THEN
2728	td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1)
2729	ENDIF
2730	tg_gdepw_0%d_value(ji,jj,ik,1) = td_risfdep%d_value(ji,jj,1,1)
2731	!gm Bug? check the gdepw_0
2732	! ... on ik
2733	tg_gdept_0%d_value(ji,jj,ik,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) &
2734	& - ( tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_0%d_value(ji,jj,ik,1) ) &
2735	& * ( tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdept_1d%d_value(1,1, ik,1) ) &
2736	& / ( tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_1d%d_value(1,1, ik,1) )
2737
2738	tg_e3t_0%d_value(ji,jj,ik ,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_0%d_value(ji,jj,ik,1)
2739	tg_e3w_0%d_value(ji,jj,ik+1,1) = tg_gdept_1d%d_value(1,1,ik+1,1) - tg_gdept_0%d_value(ji,jj,ik,1)
2740
2741	IF( ik + 1 == tg_mbathy%d_value(ji,jj,1,1) )THEN ! ice shelf point only (2 cell water column)
2742	tg_e3w_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value(ji,jj,ik+1,1) - tg_gdept_0%d_value(ji,jj,ik,1)
2743	ENDIF
2744	! ... on ik / ik-1
2745	tg_e3w_0%d_value (ji,jj,ik ,1) = tg_e3t_0%d_value (ji,jj,ik,1)
2746	tg_e3t_0%d_value (ji,jj,ik-1,1) = tg_gdepw_0%d_value(ji,jj,ik,1) - tg_gdepw_1d%d_value(1,1,ik-1,1)
2747	! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
2748	tg_gdept_0%d_value(ji,jj,ik-1,1) = tg_gdept_1d%d_value(1,1,ik-1,1)
2749
2750	ENDIF
2751	END DO
2752	END DO
2753
2754	it = 0
2755	DO jj = 1, jpj
2756	DO ji = 1, jpi
2757	ik = tg_misfdep%d_value(ji,jj,1,1)
2758	IF( ik > 1 ) THEN ! ice shelf point only
2759	tg_e3tp%d_value(ji,jj,1,1) = tg_e3t_0%d_value(ji,jj,ik ,1)
2760	tg_e3wp%d_value(ji,jj,1,1) = tg_e3w_0%d_value(ji,jj,ik+1,1)
2761	! test
2762	zdiff= tg_gdept_0%d_value(ji,jj,ik,1) - tg_gdepw_0%d_value(ji,jj,ik,1)
2763	IF( zdiff <= 0. ) THEN
2764	it = it + 1
2765	CALL logger_info(' it = '//TRIM(fct_str(it))//&
2766	& ' ik = '//TRIM(fct_str(ik))//&
2767	& ' (i,j) = '//trim(fct_str(ji))//' '//TRIM(fct_str(jj)))
2768
2769	CALL logger_info(' risfdep = '//TRIM(fct_str(td_risfdep%d_value(ji,jj,1,1))))
2770	CALL logger_info(' gdept = '//TRIM(fct_str(tg_gdept_0%d_value(ji,jj,ik,1)))//&
2771	& ' gdepw = '//TRIM(fct_str(tg_gdepw_0%d_value(ji,jj,ik+1,1)))//&
2772	& ' zdiff = '//TRIM(fct_str(zdiff)))
2773	CALL logger_info(' e3tp = '//TRIM(fct_str(tg_e3t_0%d_value(ji,jj,ik ,1)))//&
2774	& ' e3wp = '//TRIM(fct_str(tg_e3w_0%d_value(ji,jj,ik+1,1))))
2775	ENDIF
2776	ENDIF
2777	END DO
2778	END DO
2779
2780	END SUBROUTINE grid_zgr__isf_fill
2781	! !-------------------------------------------------------------------
2782	! !> @brief This subroutine define e3t, u, v, w for ISF case
2783	! !>
2784	! !> @details
2785	! !>
2786	! !> @author J.Paul
2787	! !> @date September, 2015 - rewrite from zgr_zps
2788	! !>
2789	! !> @param[in] jpi
2790	! !> @param[in] jpj
2791	! !> @param[in] td_risfdep
2792	! !-------------------------------------------------------------------
2793	! SUBROUTINE grid_zgr__isf_fill_e3x( jpi,jpj, &
2794	! & td_risfdep)
2795	! IMPLICIT NONE
2796	! ! Argument
2797	! INTEGER(i4), INTENT(IN ) :: jpi
2798	! INTEGER(i4), INTENT(IN ) :: jpj
2799	! TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
2800	!
2801	! ! local variable
2802	! REAL(dp) :: zdiff ! temporary scalar
2803	!
2804	! ! loop indices
2805	! INTEGER(i4) :: ji
2806	! INTEGER(i4) :: jj
2807	! INTEGER(i4) :: it
2808	! INTEGER(i4) :: ik
2809	! !----------------------------------------------------------------
2810	!
2811	! ! (ISF) Definition of e3t, u, v, w for ISF case
2812	! DO jj = 1, jpj
2813	! DO ji = 1, jpi
2814	! ik = tg_misfdep%d_value(ji,jj,1,1)
2815	!
2816	! IF( ik > 1 ) THEN
2817	! ! ice shelf point only
2818	! IF( td_risfdep%d_value(ji,jj,1,1) < tg_gdepw_1d%d_value(1,1,ik,1) )THEN
2819	! td_risfdep%d_value(ji,jj,1,1) = tg_gdepw_1d%d_value(1,1,ik,1)
2820	! ENDIF
2821	! tg_gdepw_0%d_value(ji,jj,ik,1) = td_risfdep%d_value(ji,jj,1,1)
2822	! !gm Bug? check the gdepw_0
2823	! ! ... on ik
2824	! tg_gdept_0%d_value(ji,jj,ik ,1) = tg_gdepw_1d%d_value(1,1,ik+1,1) - &
2825	! & (tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_0%d_value (ji,jj,ik,1)) * &
2826	! & (tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdept_1d%d_value( 1, 1,ik,1)) / &
2827	! & (tg_gdepw_1d%d_value(1,1,ik+1,1) - tg_gdepw_1d%d_value( 1, 1,ik,1))
2828	!
2829	! tg_e3t_0%d_value (ji,jj,ik ,1) = tg_gdepw_1d%d_value( 1, 1,ik+1,1) - &
2830	! & tg_gdepw_0%d_value (ji,jj,ik ,1)
2831	!
2832	! tg_e3w_0%d_value (ji,jj,ik+1,1) = tg_gdept_1d%d_value( 1, 1,ik+1,1) - &
2833	! & tg_gdept_0%d_value (ji,jj,ik ,1)
2834	!
2835	! IF( ik + 1 == tg_mbathy%d_value(ji,jj,1,1) ) THEN
2836	!
2837	! ! ice shelf point only (2 cell water column)
2838	! tg_e3w_0%d_value(ji,jj,ik+1,1) = tg_gdept_0%d_value(ji,jj,ik+1,1) - &
2839	! & tg_gdept_0%d_value(ji,jj,ik ,1)
2840	!
2841	! ENDIF
2842	! ! ... on ik / ik-1
2843	! tg_e3w_0%d_value (ji,jj,ik ,1) = 2._dp * (tg_gdept_0%d_value(ji,jj,ik,1) - &
2844	! & tg_gdepw_0%d_value(ji,jj,ik,1))
2845	!
2846	! tg_e3t_0%d_value (ji,jj,ik-1,1) = tg_gdepw_0%d_value (ji,jj,ik ,1) - &
2847	! & tg_gdepw_1d%d_value( 1, 1,ik-1,1)
2848	!
2849	! ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code
2850	! tg_gdept_0%d_value(ji,jj,ik-1,1) = tg_gdept_1d%d_value(1,1,ik-1,1)
2851	! ENDIF
2852	!
2853	! END DO
2854	! END DO
2855	!
2856	! it = 0
2857	! DO jj = 1, jpj
2858	! DO ji = 1, jpi
2859	! ik = tg_misfdep%d_value(ji,jj,1,1)
2860	! IF( ik > 1 ) THEN ! ice shelf point only
2861	! tg_e3tp%d_value(ji,jj,1,1) = tg_e3t_0%d_value(ji,jj,ik ,1)
2862	! tg_e3wp%d_value(ji,jj,1,1) = tg_e3w_0%d_value(ji,jj,ik+1,1)
2863	! ! test
2864	! zdiff= tg_gdept_0%d_value(ji,jj,ik,1) - &
2865	! & tg_gdepw_0%d_value(ji,jj,ik,1)
2866	!
2867	! IF( zdiff <= 0. ) THEN
2868	! it = it + 1
2869	! CALL logger_info(' it = '//TRIM(fct_str(it))//&
2870	! & ' ik = '//TRIM(fct_str(ik))//&
2871	! & ' (i,j) =('//TRIM(fct_str(ji))//','//TRIM(fct_str(jj))//')')
2872	! CALL logger_info(' risfdep = '//TRIM(fct_str(td_risfdep%d_value(ji,jj,1,1))) )
2873	! CALL logger_info(' gdept = '//TRIM(fct_str(tg_gdept_0%d_value(ji,jj,ik ,1)))//&
2874	! & ' gdepw = '//TRIM(fct_str(tg_gdepw_0%d_value(ji,jj,ik+1,1)))//&
2875	! & ' zdiff = '//TRIM(fct_str(zdiff)) )
2876	! CALL logger_info(' e3tp = '//TRIM(fct_str( tg_e3tp%d_value(ji,jj,1,1)))//&
2877	! & ' e3wp = '//TRIM(fct_str( tg_e3wp%d_value(ji,jj,1,1))) )
2878	! ENDIF
2879	! ENDIF
2880	! END DO
2881	! END DO
2882	! ! END (ISF)
2883	!
2884	! END SUBROUTINE grid_zgr__isf_fill_e3x
2885	!-------------------------------------------------------------------
2886	!> @brief This subroutine define e3uw
2887	!> (adapted for 2 cells in the water column) for ISF case
2888	!>
2889	!> @details
2890	!>
2891	!> @author J.Paul
2892	!> @date September, 2015 - rewrite from zgr_zps
2893	!>
2894	!> @param[in] jpi
2895	!> @param[in] jpj
2896	!-------------------------------------------------------------------
2897	SUBROUTINE grid_zgr__isf_fill_e3uw(jpi,jpj)
2898	IMPLICIT NONE
2899	! Argument
2900	INTEGER(i4) , INTENT(IN ) :: jpi
2901	INTEGER(i4) , INTENT(IN ) :: jpj
2902
2903	! local variable
2904	INTEGER(i4) :: ikb, ikt
2905
2906	! loop indices
2907	INTEGER(i4) :: ji
2908	INTEGER(i4) :: jj
2909	!----------------------------------------------------------------
2910
2911	DO jj = 2, jpj - 1
2912	DO ji = 2, jpi - 1
2913
2914	ikb = MAX(tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji+1,jj,1,1))
2915	ikt = MAX(tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji+1,jj,1,1))
2916	IF( ikb == ikt+1 )THEN
2917	tg_e3uw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb ,1), tg_gdept_0%d_value(ji+1,jj ,ikb ,1) ) - &
2918	& MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji+1,jj ,ikb-1,1) )
2919	ENDIF
2920
2921	ikb = MAX( tg_mbathy%d_value (ji,jj,1,1), tg_mbathy%d_value (ji,jj+1,1,1))
2922	ikt = MAX( tg_misfdep%d_value(ji,jj,1,1), tg_misfdep%d_value(ji,jj+1,1,1))
2923	IF( ikb == ikt+1 )THEN
2924	tg_e3vw_0%d_value(ji,jj,ikb,1) = MIN( tg_gdept_0%d_value(ji,jj,ikb ,1), tg_gdept_0%d_value(ji ,jj+1,ikb ,1) ) - &
2925	& MAX( tg_gdept_0%d_value(ji,jj,ikb-1,1), tg_gdept_0%d_value(ji ,jj+1,ikb-1,1) )
2926	ENDIF
2927	END DO
2928	END DO
2929
2930	END SUBROUTINE grid_zgr__isf_fill_e3uw
2931	! !-------------------------------------------------------------------
2932	! !> @brief This subroutine compute gdep3w_0 (vertical sum of e3w)
2933	! !>
2934	! !> @details
2935	! !>
2936	! !> @author J.Paul
2937	! !> @date September, 2015 - rewrite from zgr_zps
2938	! !>
2939	! !> @param[in] jpi
2940	! !> @param[in] jpj
2941	! !> @param[in] jpk
2942	! !> @param[in] td_risfdep
2943	! !-------------------------------------------------------------------
2944	! SUBROUTINE grid_zgr__isf_fill_gdep3w_0( jpi,jpj,jpk,td_risfdep )
2945	! IMPLICIT NONE
2946	! ! Argument
2947	! INTEGER(i4), INTENT(IN ) :: jpi
2948	! INTEGER(i4), INTENT(IN ) :: jpj
2949	! INTEGER(i4), INTENT(IN ) :: jpk
2950	! TYPE(TVAR) , INTENT(INOUT) :: td_risfdep
2951	!
2952	! ! local variable
2953	! INTEGER(i4) :: ik
2954	!
2955	! ! loop indices
2956	! INTEGER(i4) :: ji
2957	! INTEGER(i4) :: jj
2958	! INTEGER(i4) :: jk
2959	! !----------------------------------------------------------------
2960	!
2961	! WHERE( tg_misfdep%d_value(:,:,:,:) == 0 ) tg_misfdep%d_value(:,:,:,:) = 1
2962	!
2963	! DO jj = 1,jpj
2964	! DO ji = 1,jpi
2965	!
2966	! tg_gdep3w_0%d_value(ji,jj,1,1) = 0.5_dp * tg_e3w_0%d_value(ji,jj,1,1)
2967	! DO jk = 2, INT(tg_misfdep%d_value(ji,jj,1,1),i4)
2968	! tg_gdep3w_0%d_value(ji,jj,jk,1) = tg_gdep3w_0%d_value(ji,jj,jk-1,1) + &
2969	! & tg_e3w_0%d_value (ji,jj,jk ,1)
2970	! END DO
2971	!
2972	! ik=tg_misfdep%d_value(ji,jj,1,1)
2973	! IF( ik >= 2 )THEN
2974	! tg_gdep3w_0%d_value(ji,jj,ik,1) = td_risfdep%d_value(ji,jj, 1,1) + &
2975	! & 0.5_dp * tg_e3w_0%d_value(ji,jj,ik,1)
2976	! ENDIF
2977	!
2978	! DO jk = ik + 1, jpk
2979	! tg_gdep3w_0%d_value(ji,jj,jk,1) = tg_gdep3w_0%d_value(ji,jj,jk-1,1) + &
2980	! & tg_e3w_0%d_value (ji,jj,jk ,1)
2981	! END DO
2982	!
2983	! END DO
2984	! END DO
2985	!
2986	! END SUBROUTINE grid_zgr__isf_fill_gdep3w_0
2987	!-------------------------------------------------------------------
2988	!> @brief This function initialise global variable needed to compute vertical
2989	!> mesh
2990	!>
2991	!> @author J.Paul
2992	!> @date September, 2015 - Initial version
2993	!>
2994	!> @param[in] jpi
2995	!> @param[in] jpj
2996	!-------------------------------------------------------------------
2997	SUBROUTINE grid_zgr_sco_init(jpi,jpj)
2998	IMPLICIT NONE
2999	! Argument
3000	INTEGER(i4), INTENT(IN) :: jpi
3001	INTEGER(i4), INTENT(IN) :: jpj
3002
3003	! local variable
3004	REAL(dp), DIMENSION(jpi,jpj) :: dl_tmp
3005
3006	! loop indices
3007	!----------------------------------------------------------------
3008
3009	dl_tmp(:,:)=dp_fill
3010
3011	tg_rx1 =var_init('rx1 ',dl_tmp(:,:), dd_fill=dp_fill, id_type=NF90_DOUBLE)
3012
3013	END SUBROUTINE grid_zgr_sco_init
3014	!-------------------------------------------------------------------
3015	!> @brief This function clean structure
3016	!>
3017	!> @author J.Paul
3018	!> @date September, 2015 - Initial version
3019	!>
3020	!-------------------------------------------------------------------
3021	SUBROUTINE grid_zgr_sco_clean()
3022	IMPLICIT NONE
3023	! Argument
3024
3025	! local variable
3026	! loop indices
3027	!----------------------------------------------------------------
3028
3029	CALL var_clean(tg_rx1 )
3030
3031	END SUBROUTINE grid_zgr_sco_clean
3032	!-------------------------------------------------------------------
3033	!> @brief This subroutine define the s-coordinate system
3034	!>
3035	!> @details
3036	!>
3037	!> ** Method : s-coordinate
3038	!> The depth of model levels is defined as the product of an
3039	!> analytical function by the local bathymetry, while the vertical
3040	!> scale factors are defined as the product of the first derivative
3041	!> of the analytical function by the bathymetry.
3042	!> (this solution save memory as depth and scale factors are not
3043	!> 3d fields)
3044	!> - Read bathymetry (in meters) at t-point and compute the
3045	!> bathymetry at u-, v-, and f-points.
3046	!> - hbatu = mi( hbatt )
3047	!> - hbatv = mj( hbatt )
3048	!> - hbatf = mi( mj( hbatt ) )
3049	!> - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
3050	!> function and its derivative given as function.
3051	!> - z_gsigt(k) = fssig (k )
3052	!> - z_gsigw(k) = fssig (k-0.5)
3053	!> - z_esigt(k) = fsdsig(k )
3054	!> - z_esigw(k) = fsdsig(k-0.5)
3055	!> Three options for stretching are give, and they can be modified
3056	!> following the users requirements. Nevertheless, the output as
3057	!> well as the way to compute the model levels and scale factors
3058	!> must be respected in order to insure second order accuracy
3059	!> schemes.
3060	!>
3061	!> The three methods for stretching available are:
3062	!>
3063	!> s_sh94 (Song and Haidvogel 1994)
3064	!> a sinh/tanh function that allows sigma and stretched sigma
3065	!>
3066	!> s_sf12 (Siddorn and Furner 2012?)
3067	!> allows the maintenance of fixed surface and or
3068	!> bottom cell resolutions (cf. geopotential coordinates)
3069	!> within an analytically derived stretched S-coordinate framework.
3070	!>
3071	!> s_tanh (Madec et al 1996)
3072	!> a cosh/tanh function that gives stretched coordinates
3073	!>
3074	!> @author J.Paul
3075	!> @date September, 2015 - rewrite from zgr_sco
3076	!> @date October, 2016
3077	!> - add wetting and drying boolean
3078	!>
3079	!-------------------------------------------------------------------
3080	SUBROUTINE grid_zgr__sco_fill( td_nam,jpi,jpj,jpk,td_bathy )
3081	IMPLICIT NONE
3082	! Argument
3083	TYPE(TNAMZ), INTENT(IN ) :: td_nam
3084	INTEGER(i4), INTENT(IN ) :: jpi
3085	INTEGER(i4), INTENT(IN ) :: jpj
3086	INTEGER(i4), INTENT(IN ) :: jpk
3087	TYPE(TVAR) , INTENT(INOUT) :: td_bathy
3088
3089	! local variable
3090	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
3091
3092	REAL(dp) :: zrmax, zrfact
3093	REAL(dp) :: ztaper
3094
3095	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hifv !< interface depth between stretching at v--f
3096	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hiff
3097	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hift !< and quasi-uniform spacing t--u points (m)
3098	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hifu
3099	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_scosrf !< ocean surface topography
3100	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_scobot !< ocean bottom topography
3101
3102	! REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hbatt
3103	! REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hbatu
3104	! REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hbatv
3105	! REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_hbatf
3106
3107	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zenv
3108	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zri
3109	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zrj
3110	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zhbat
3111	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpi1
3112	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpi2
3113	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpj1
3114	REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ztmpj2
3115
3116	! loop indices
3117	INTEGER(i4) :: ji
3118	INTEGER(i4) :: jj
3119	INTEGER(i4) :: jk
3120	INTEGER(i4) :: jl
3121	!----------------------------------------------------------------
3122
3123	CALL logger_info('GRID ZGR SCO: s-coordinate or hybrid z-s-coordinate')
3124	CALL logger_info('~~~~~~~~~~~')
3125	CALL logger_info(' Namelist namzgr_sco')
3126	CALL logger_info(' stretching coeffs ')
3127	CALL logger_info(' maximum depth of s-bottom surface (>0) dn_sbot_max = '//TRIM(fct_str(td_nam%d_sbot_max)))
3128	CALL logger_info(' minimum depth of s-bottom surface (>0) dn_sbot_min = '//TRIM(fct_str(td_nam%d_sbot_min)))
3129	CALL logger_info(' Critical depth dn_hc = '//TRIM(fct_str(td_nam%d_hc)))
3130	CALL logger_info(' maximum cut-off r-value allowed dn_rmax = '//TRIM(fct_str(td_nam%d_rmax)))
3131	CALL logger_info(' Song and Haidvogel 1994 stretching ln_s_sh94 = '//TRIM(fct_str(td_nam%l_s_sh94)))
3132	CALL logger_info(' Song and Haidvogel 1994 stretching coefficients')
3133	CALL logger_info(' surface control parameter (0<=rn_theta<=20) dn_theta = '//TRIM(fct_str(td_nam%d_theta)))
3134	CALL logger_info(' bottom control parameter (0<=rn_thetb<= 1) dn_thetb = '//TRIM(fct_str(td_nam%d_thetb)))
3135	CALL logger_info(' stretching parameter (song and haidvogel) dn_bb = '//TRIM(fct_str(td_nam%d_bb)))
3136	CALL logger_info(' Siddorn and Furner 2012 stretching ln_s_sf12 = '//TRIM(fct_str(td_nam%l_s_sf12)))
3137	CALL logger_info(' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = '//TRIM(fct_str(td_nam%l_sigcrit)))
3138	CALL logger_info(' Siddorn and Furner 2012 stretching coefficients')
3139	CALL logger_info(' stretchin parameter ( >1 surface; <1 bottom) dn_alpha = '//TRIM(fct_str(td_nam%d_alpha)))
3140	CALL logger_info(' e-fold length scale for transition region dn_efold = '//TRIM(fct_str(td_nam%d_efold)))
3141	CALL logger_info(' Surface cell depth (Zs) (m) dn_zs = '//TRIM(fct_str(td_nam%d_zs)))
3142	CALL logger_info(' Bathymetry multiplier for Zb dn_zb_a = '//TRIM(fct_str(td_nam%d_zb_a)))
3143	CALL logger_info(' Offset for Zb dn_zb_b = '//TRIM(fct_str(td_nam%d_zb_b)))
3144	CALL logger_info(' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b')
3145
3146	ALLOCATE(dl_hifv(jpi,jpj))
3147	ALLOCATE(dl_hiff(jpi,jpj))
3148	ALLOCATE(dl_hift(jpi,jpj))
3149	ALLOCATE(dl_hifu(jpi,jpj))
3150
3151	! set the minimum depth for the s-coordinate
3152	dl_hift(:,:) = td_nam%d_sbot_min
3153	dl_hifu(:,:) = td_nam%d_sbot_min
3154	dl_hifv(:,:) = td_nam%d_sbot_min
3155	dl_hiff(:,:) = td_nam%d_sbot_min
3156
3157	! set maximum ocean depth
3158	td_bathy%d_value(:,:,1,1) = MIN( td_nam%d_sbot_max, td_bathy%d_value(:,:,1,1) )
3159	IF( .NOT. td_nam%l_wd )THEN
3160	DO jj = 1, jpj
3161	DO ji = 1, jpi
3162	IF( td_bathy%d_value(ji,jj,1,1) > 0._dp )THEN
3163	td_bathy%d_value(ji,jj,1,1) = MAX(td_nam%d_sbot_min, td_bathy%d_value(ji,jj,1,1))
3164	ENDIF
3165	END DO
3166	END DO
3167	ENDIF
3168
3169	! =============================
3170	! Define the envelop bathymetry (hbatt)
3171	! =============================
3172	! use r-value to create hybrid coordinates
3173	ALLOCATE(zenv(jpi,jpj))
3174	zenv(:,:) = td_bathy%d_value(:,:,1,1)
3175
3176	IF( .NOT. td_nam%l_wd )THEN
3177	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
3178	DO jj = 1, jpj
3179	DO ji = 1, jpi
3180	IF( td_bathy%d_value(ji,jj,1,1) == 0._dp )THEN
3181	iip1 = MIN( ji+1, jpi )
3182	ijp1 = MIN( jj+1, jpj )
3183	iim1 = MAX( ji-1, 1 )
3184	ijm1 = MAX( jj-1, 1 )
3185	IF( ( td_bathy%d_value(iim1,ijm1,1,1) + &
3186	& td_bathy%d_value(ji ,ijp1,1,1) + &
3187	& td_bathy%d_value(iip1,ijp1,1,1) + &
3188	& td_bathy%d_value(iim1,jj ,1,1) + &
3189	& td_bathy%d_value(iip1,jj ,1,1) + &
3190	& td_bathy%d_value(iim1,ijm1,1,1) + &
3191	& td_bathy%d_value(ji ,ijm1,1,1) + &
3192	& td_bathy%d_value(iip1,ijp1,1,1) ) > 0._dp )THEN
3193	zenv(ji,jj) = td_nam%d_sbot_min
3194	ENDIF
3195	ENDIF
3196	END DO
3197	END DO
3198	ENDIF
3199
3200	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
3201	! this is only in mpp case, so here just do nothing
3202	!! CALL lbc_lnk( zenv(:,:), 'T', td_nam%i_perio, 1._dp, 'no0' )
3203
3204	! smooth the bathymetry (if required)
3205	ALLOCATE(dl_scosrf(jpi,jpj))
3206	ALLOCATE(dl_scobot(jpi,jpj))
3207	dl_scosrf(:,:) = 0._dp ! ocean surface depth (here zero: no under ice-shelf sea)
3208	dl_scobot(:,:) = td_bathy%d_value(:,:,1,1) ! ocean bottom depth
3209
3210	jl = 0
3211	zrmax = 1._dp
3212
3213	! set scaling factor used in reducing vertical gradients
3214	zrfact = ( 1._dp - td_nam%d_rmax ) / ( 1._dp + td_nam%d_rmax )
3215
3216	! initialise temporary envelope depth arrays
3217	ALLOCATE(ztmpi1(jpi,jpj))
3218	ALLOCATE(ztmpi2(jpi,jpj))
3219	ALLOCATE(ztmpj1(jpi,jpj))
3220	ALLOCATE(ztmpj2(jpi,jpj))
3221
3222	ztmpi1(:,:) = zenv(:,:)
3223	ztmpi2(:,:) = zenv(:,:)
3224	ztmpj1(:,:) = zenv(:,:)
3225	ztmpj2(:,:) = zenv(:,:)
3226
3227	! initialise temporary r-value arrays
3228	ALLOCATE(zri(jpi,jpj))
3229	ALLOCATE(zrj(jpi,jpj))
3230	zri(:,:) = 1._dp
3231	zrj(:,:) = 1._dp
3232
3233	! Iterative loop !
3234	! ================ !
3235	DO WHILE( jl <= 10000 .AND. ( zrmax - td_nam%d_rmax ) > 1.e-8_dp )
3236
3237	jl = jl + 1
3238	zrmax = 0._dp
3239	! we set zrmax from previous r-values (zri and zrj) first
3240	! if set after current r-value calculation (as previously)
3241	! we could exit DO WHILE prematurely before checking r-value
3242	! of current zenv
3243	DO jj = 1, jpj
3244	DO ji = 1, jpi
3245	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
3246	END DO
3247	END DO
3248	zri(:,:) = 0._dp
3249	zrj(:,:) = 0._dp
3250	DO jj = 1, jpj
3251	DO ji = 1, jpi
3252	iip1 = MIN( ji+1, jpi ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
3253	ijp1 = MIN( jj+1, jpj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
3254	IF( (zenv(ji ,jj) > 0._dp) .AND. &
3255	& (zenv(iip1,jj) > 0._dp) )THEN
3256	zri(ji,jj) = ( zenv(iip1, jj) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
3257	END IF
3258	IF( (zenv(ji,jj ) > 0._dp) .AND. &
3259	& (zenv(ji,ijp1) > 0._dp) )THEN
3260	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
3261	END IF
3262	IF( zri(ji,jj) > td_nam%d_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
3263	IF( zri(ji,jj) < -td_nam%d_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
3264	IF( zrj(ji,jj) > td_nam%d_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
3265	IF( zrj(ji,jj) < -td_nam%d_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
3266	END DO
3267	END DO
3268	!!
3269	!
3270	CALL logger_info('zgr_sco : iter= '//TRIM(fct_str(jl))//&
3271	& ' rmax= '//TRIM(fct_str(zrmax)) )
3272	!
3273	DO jj = 1, jpj
3274	DO ji = 1, jpi
3275	zenv(ji,jj) = MAX( zenv (ji,jj), &
3276	& ztmpi1(ji,jj), &
3277	& ztmpi2(ji,jj), &
3278	& ztmpj1(ji,jj), &
3279	& ztmpj2(ji,jj) )
3280	END DO
3281	END DO
3282	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
3283	! this is only in mpp case, so here just do nothing
3284	!!CALL lbc_lnk( zenv, 'T', td_nam%i_perio, 1._dp, 'no0' )
3285	ENDDO
3286	! End loop !
3287	! ================ !
3288
3289	DEALLOCATE(zri)
3290	DEALLOCATE(zrj)
3291
3292	DEALLOCATE(ztmpi1)
3293	DEALLOCATE(ztmpi2)
3294	DEALLOCATE(ztmpj1)
3295	DEALLOCATE(ztmpj2)
3296
3297	DO jj = 1, jpj
3298	DO ji = 1, jpi
3299	! set all points to avoid undefined scale value warnings
3300	zenv(ji,jj) = MAX( zenv(ji,jj), td_nam%d_sbot_min )
3301	ENDDO
3302	ENDDO
3303	!
3304	! Envelope bathymetry saved in hbatt
3305	!
3306	! ALLOCATE(dl_hbatt(jpi,jpj))
3307	! ALLOCATE(dl_hbatu(jpi,jpj))
3308	! ALLOCATE(dl_hbatv(jpi,jpj))
3309	! ALLOCATE(dl_hbatf(jpi,jpj))
3310
3311	tg_hbatt%d_value(:,:,1,1) = zenv(:,:)
3312	IF( MINVAL( tg_gphit%d_value(:,:,1,1) ) * &
3313	& MAXVAL( tg_gphit%d_value(:,:,1,1) ) <= 0._dp ) THEN
3314	CALL logger_warn( ' s-coordinates are tapered in vicinity of the Equator' )
3315	DO jj = 1, jpj
3316	DO ji = 1, jpi
3317	ztaper = EXP( -(tg_gphit%d_value(ji,jj,1,1)/8._dp)**2._dp )
3318	tg_hbatt%d_value(ji,jj,1,1) = td_nam%d_sbot_max * ztaper &
3319	& + tg_hbatt%d_value(ji,jj,1,1) * (1._dp - ztaper)
3320	END DO
3321	END DO
3322	ENDIF
3323
3324	DEALLOCATE(zenv)
3325
3326	CALL logger_info(' bathy MAX '//TRIM(fct_str(MAXVAL( td_bathy%d_value(:,:,1,1) )))//&
3327	& ' MIN '//TRIM(fct_str(MINVAL( td_bathy%d_value(:,:,1,1) ))) )
3328	CALL logger_info(' hbatt MAX '//TRIM(fct_str(MAXVAL( tg_hbatt%d_value(:,:,1,1) )))//&
3329	& ' MIN '//TRIM(fct_str(MINVAL( tg_hbatt%d_value(:,:,1,1) ))) )
3330	!
3331	! hbatu, hbatv, hbatf fields
3332	!
3333	tg_hbatu%d_value(:,:,1,1) = td_nam%d_sbot_min
3334	tg_hbatv%d_value(:,:,1,1) = td_nam%d_sbot_min
3335	tg_hbatf%d_value(:,:,1,1) = td_nam%d_sbot_min
3336
3337	DO jj = 1, jpj-1
3338	DO ji = 1, jpi-1 ! NO vector opt.
3339	tg_hbatu%d_value(ji,jj,1,1) = 0.50_dp * ( tg_hbatt%d_value(ji ,jj ,1,1) + &
3340	& tg_hbatt%d_value(ji+1,jj ,1,1) )
3341	tg_hbatv%d_value(ji,jj,1,1) = 0.50_dp * ( tg_hbatt%d_value(ji ,jj ,1,1) + &
3342	& tg_hbatt%d_value(ji ,jj+1,1,1) )
3343	tg_hbatf%d_value(ji,jj,1,1) = 0.25_dp * ( tg_hbatt%d_value(ji ,jj ,1,1) + &
3344	& tg_hbatt%d_value(ji ,jj+1,1,1) + &
3345	& tg_hbatt%d_value(ji+1,jj ,1,1) + &
3346	& tg_hbatt%d_value(ji+1,jj+1,1,1) )
3347	ENDDO
3348	ENDDO
3349
3350	IF( td_nam%l_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points
3351	DO jj = 1, jpj
3352	DO ji = 1, jpi
3353	IF( ABS(tg_hbatt%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
3354	tg_hbatt%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatt%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
3355	ENDIF
3356	IF( ABS(tg_hbatu%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
3357	tg_hbatu%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatu%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
3358	ENDIF
3359	IF( ABS(tg_hbatv%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
3360	tg_hbatv%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatv%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
3361	ENDIF
3362	IF( ABS(tg_hbatf%d_value(ji,jj,1,1)) < td_nam%d_wdmin1 )THEN
3363	tg_hbatf%d_value(ji,jj,1,1) = SIGN(1._dp, tg_hbatf%d_value(ji,jj,1,1)) * td_nam%d_wdmin1
3364	ENDIF
3365	END DO
3366	END DO
3367	ENDIF
3368
3369	! Apply lateral boundary condition
3370	ALLOCATE(zhbat(jpi,jpj))
3371
3372	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
3373	zhbat(:,:) = tg_hbatu%d_value(:,:,1,1)
3374	CALL lbc_lnk( tg_hbatu%d_value(:,:,1,1), 'U', td_nam%i_perio, 1._dp )
3375	DO jj = 1, jpj
3376	DO ji = 1, jpi
3377	IF( tg_hbatu%d_value(ji,jj,1,1) == 0._dp )THEN
3378	!No worries about the following line when l_wd == .true.
3379	IF( zhbat(ji,jj) == 0._dp ) tg_hbatu%d_value(ji,jj,1,1) = td_nam%d_sbot_min
3380	IF( zhbat(ji,jj) /= 0._dp ) tg_hbatu%d_value(ji,jj,1,1) = zhbat(ji,jj)
3381	ENDIF
3382	ENDDO
3383	ENDDO
3384
3385	zhbat(:,:) = tg_hbatv%d_value(:,:,1,1)
3386	CALL lbc_lnk( tg_hbatv%d_value(:,:,1,1), 'V', td_nam%i_perio, 1._dp )
3387	DO jj = 1, jpj
3388	DO ji = 1, jpi
3389	IF( tg_hbatv%d_value(ji,jj,1,1) == 0._dp ) THEN
3390	IF( zhbat(ji,jj) == 0._dp ) tg_hbatv%d_value(ji,jj,1,1) = td_nam%d_sbot_min
3391	IF( zhbat(ji,jj) /= 0._dp ) tg_hbatv%d_value(ji,jj,1,1) = zhbat(ji,jj)
3392	ENDIF
3393	ENDDO
3394	ENDDO
3395
3396	zhbat(:,:) = tg_hbatf%d_value(:,:,1,1)
3397	CALL lbc_lnk( tg_hbatf%d_value(:,:,1,1), 'F', td_nam%i_perio, 1._dp )
3398	DO jj = 1, jpj
3399	DO ji = 1, jpi
3400	IF( tg_hbatf%d_value(ji,jj,1,1) == 0._dp ) THEN
3401	IF( zhbat(ji,jj) == 0._dp ) tg_hbatf%d_value(ji,jj,1,1) = td_nam%d_sbot_min
3402	IF( zhbat(ji,jj) /= 0._dp ) tg_hbatf%d_value(ji,jj,1,1) = zhbat(ji,jj)
3403	ENDIF
3404	ENDDO
3405	ENDDO
3406
3407	DEALLOCATE(zhbat)
3408
3409	!!bug: key_helsinki a verifer
3410	IF( .NOT.td_nam%l_wd )THEN
3411	dl_hift(:,:) = MIN( dl_hift(:,:), tg_hbatt%d_value(:,:,1,1) )
3412	dl_hifu(:,:) = MIN( dl_hifu(:,:), tg_hbatu%d_value(:,:,1,1) )
3413	dl_hifv(:,:) = MIN( dl_hifv(:,:), tg_hbatv%d_value(:,:,1,1) )
3414	dl_hiff(:,:) = MIN( dl_hiff(:,:), tg_hbatf%d_value(:,:,1,1) )
3415	ENDIF
3416
3417	CALL logger_info(' MAX val hif t '//TRIM(fct_str(MAXVAL( dl_hift(:,:) )))//&
3418	& ' f '//TRIM(fct_str(MAXVAL( dl_hiff(:,:) )))//&
3419	& ' u '//TRIM(fct_str(MAXVAL( dl_hifu(:,:) )))//&
3420	& ' v '//TRIM(fct_str(MAXVAL( dl_hifv(:,:) ))) )
3421	CALL logger_info(' MIN val hif t '//TRIM(fct_str(MINVAL( dl_hift(:,:) )))//&
3422	& ' f '//TRIM(fct_str(MINVAL( dl_hiff(:,:) )))//&
3423	& ' u '//TRIM(fct_str(MINVAL( dl_hifu(:,:) )))//&
3424	& ' v '//TRIM(fct_str(MINVAL( dl_hifv(:,:) ))) )
3425	CALL logger_info(' MAX val hbat t '//TRIM(fct_str(MAXVAL( tg_hbatt%d_value(:,:,1,1) )))//&
3426	& ' f '//TRIM(fct_str(MAXVAL( tg_hbatf%d_value(:,:,1,1) )))//&
3427	& ' u '//TRIM(fct_str(MAXVAL( tg_hbatu%d_value(:,:,1,1) )))//&
3428	& ' v '//TRIM(fct_str(MAXVAL( tg_hbatv%d_value(:,:,1,1) ))) )
3429	CALL logger_info(' MIN val hbat t '//TRIM(fct_str(MINVAL( tg_hbatt%d_value(:,:,1,1) )))//&
3430	& ' f '//TRIM(fct_str(MINVAL( tg_hbatf%d_value(:,:,1,1) )))//&
3431	& ' u '//TRIM(fct_str(MINVAL( tg_hbatu%d_value(:,:,1,1) )))//&
3432	& ' v '//TRIM(fct_str(MINVAL( tg_hbatv%d_value(:,:,1,1) ))) )
3433	!! helsinki
3434
3435	! =======================
3436	! s-ccordinate fields (gdep., e3.)
3437	! =======================
3438	!
3439	! non-dimensional "sigma" for model level depth at w- and t-levels
3440
3441	!========================================================================
3442	! Song and Haidvogel 1994 (ln_s_sh94=T)
3443	! Siddorn and Furner 2012 (ln_sf12=T)
3444	! or tanh function (both false)
3445	!========================================================================
3446	IF( td_nam%l_s_sh94 ) THEN
3447	CALL grid_zgr__sco_s_sh94( td_nam,jpi,jpj,jpk, &
3448	& dl_scosrf )
3449	ELSEIF( td_nam%l_s_sf12 ) THEN
3450	CALL grid_zgr__sco_s_sf12( td_nam,jpi,jpj,jpk, &
3451	& dl_scosrf )
3452	ELSE
3453	CALL grid_zgr__sco_s_tanh( td_nam,jpi,jpj,jpk, &
3454	& dl_scosrf, &
3455	& dl_hift, dl_hifu, dl_hifv, dl_hiff )
3456	ENDIF
3457
3458	DEALLOCATE(dl_scosrf)
3459
3460	DEALLOCATE(dl_hifv)
3461	DEALLOCATE(dl_hiff)
3462	DEALLOCATE(dl_hift)
3463	DEALLOCATE(dl_hifu)
3464
3465	CALL lbc_lnk( tg_e3t_0%d_value(:,:,:,1) , 'T', td_nam%i_perio, 1._dp )
3466	CALL lbc_lnk( tg_e3u_0%d_value(:,:,:,1) , 'U', td_nam%i_perio, 1._dp )
3467	CALL lbc_lnk( tg_e3v_0%d_value(:,:,:,1) , 'V', td_nam%i_perio, 1._dp )
3468	CALL lbc_lnk( tg_e3f_0%d_value(:,:,:,1) , 'F', td_nam%i_perio, 1._dp )
3469	CALL lbc_lnk( tg_e3w_0%d_value(:,:,:,1) , 'W', td_nam%i_perio, 1._dp )
3470	CALL lbc_lnk( tg_e3uw_0%d_value(:,:,:,1), 'U', td_nam%i_perio, 1._dp )
3471	CALL lbc_lnk( tg_e3vw_0%d_value(:,:,:,1), 'V', td_nam%i_perio, 1._dp )
3472
3473	IF( .NOT. td_nam%l_wd ) THEN
3474	WHERE( tg_e3t_0%d_value(:,:,:,1) == 0_dp ) tg_e3t_0%d_value(:,:,:,1) = 1._dp
3475	WHERE( tg_e3u_0%d_value(:,:,:,1) == 0_dp ) tg_e3u_0%d_value(:,:,:,1) = 1._dp
3476	WHERE( tg_e3v_0%d_value(:,:,:,1) == 0_dp ) tg_e3v_0%d_value(:,:,:,1) = 1._dp
3477	WHERE( tg_e3f_0%d_value(:,:,:,1) == 0_dp ) tg_e3f_0%d_value(:,:,:,1) = 1._dp
3478	WHERE( tg_e3w_0%d_value(:,:,:,1) == 0_dp ) tg_e3w_0%d_value(:,:,:,1) = 1._dp
3479	WHERE( tg_e3uw_0%d_value(:,:,:,1)== 0_dp ) tg_e3uw_0%d_value(:,:,:,1)= 1._dp
3480	WHERE( tg_e3vw_0%d_value(:,:,:,1)== 0_dp ) tg_e3vw_0%d_value(:,:,:,1)= 1._dp
3481	ENDIF
3482
3483	! HYBRID :
3484	DO jj = 1, jpj
3485	DO ji = 1, jpi
3486	DO jk = 1, jpk-1
3487	IF( dl_scobot(ji,jj) >= tg_gdept_0%d_value(ji,jj,jk,1) )THEN
3488	tg_mbathy%d_value(ji,jj,1,1) = REAL(MAX( 2, jk ),dp)
3489	ENDIF
3490	END DO
3491
3492	IF( td_nam%l_wd ) THEN
3493	IF( dl_scobot(ji,jj) <= -(td_nam%d_wdld - td_nam%d_wdmin2) )THEN
3494
3495	tg_mbathy%d_value(ji,jj,1,1) = 0._dp
3496
3497	ELSEIF( dl_scobot(ji,jj) <= td_nam%d_wdmin1 )THEN
3498
3499	tg_mbathy%d_value(ji,jj,1,1) = 2._dp
3500
3501	ENDIF
3502	ELSE
3503	IF( dl_scobot(ji,jj) == 0._dp )THEN
3504	tg_mbathy%d_value(ji,jj,1,1) = 0._dp
3505	ENDIF
3506	ENDIF
3507
3508	ENDDO
3509	ENDDO
3510
3511	DEALLOCATE(dl_scobot)
3512
3513	CALL logger_info(' MIN val mbathy MIN '//TRIM(fct_str(MINVAL( tg_mbathy%d_value(:,:,1,1) )))//&
3514	& ' MAX '//TRIM(fct_str(MAXVAL( tg_mbathy%d_value(:,:,1,1) ))) )
3515	CALL logger_info(' MIN val depth t '//TRIM(fct_str(MINVAL( tg_gdept_0%d_value(:,:,:,1) )))//&
3516	& ' w '//TRIM(fct_str(MINVAL( tg_gdepw_0%d_value(:,:,:,1) ))) )!//&
3517	!& '3w '//TRIM(fct_str(MINVAL( tg_gdep3w_0%d_value(:,:,:,1)))) )
3518
3519	CALL logger_info(' MIN val e3 t '//TRIM(fct_str(MINVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
3520	& ' f '//TRIM(fct_str(MINVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
3521	& ' u '//TRIM(fct_str(MINVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
3522	& ' v '//TRIM(fct_str(MINVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
3523	& ' uw '//TRIM(fct_str(MINVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
3524	& ' vw '//TRIM(fct_str(MINVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
3525	& ' w '//TRIM(fct_str(MINVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
3526
3527	CALL logger_info(' MAX val depth t '//TRIM(fct_str(MAXVAL( tg_gdept_0%d_value(:,:,:,1) )))//&
3528	& ' w '//TRIM(fct_str(MAXVAL( tg_gdepw_0%d_value(:,:,:,1) ))) )!//&
3529	!& '3w '//TRIM(fct_str(MAXVAL( tg_gdep3w_0%d_value(:,:,:,1) ))) )
3530
3531	CALL logger_info(' MAX val e3 t '//TRIM(fct_str(MAXVAL( tg_e3t_0%d_value(:,:,:,1) )))//&
3532	& ' f '//TRIM(fct_str(MAXVAL( tg_e3f_0%d_value(:,:,:,1) )))//&
3533	& ' u '//TRIM(fct_str(MAXVAL( tg_e3u_0%d_value(:,:,:,1) )))//&
3534	& ' v '//TRIM(fct_str(MAXVAL( tg_e3v_0%d_value(:,:,:,1) )))//&
3535	& ' uw '//TRIM(fct_str(MAXVAL( tg_e3uw_0%d_value(:,:,:,1) )))//&
3536	& ' vw '//TRIM(fct_str(MAXVAL( tg_e3vw_0%d_value(:,:,:,1) )))//&
3537	& ' w '//TRIM(fct_str(MAXVAL( tg_e3w_0%d_value(:,:,:,1) ))) )
3538
3539	!================================================================================
3540	! check the coordinate makes sense
3541	!================================================================================
3542	DO ji = 1, jpi
3543	DO jj = 1, jpj
3544
3545	IF( tg_hbatt%d_value(ji,jj,1,1) > 0._dp )THEN
3546	DO jk = 1, INT(tg_mbathy%d_value(ji,jj,1,1),i4)
3547	! check coordinate is monotonically increasing
3548	IF( tg_e3w_0%d_value(ji,jj,jk,1) <= 0._dp .OR. &
3549	& tg_e3t_0%d_value(ji,jj,jk,1) <= 0._dp )THEN
3550	CALL logger_fatal(' GRID ZGR SCO: e3w or e3t =< 0 at point ('//&
3551	& TRIM(fct_str(ji))//","//&
3552	& TRIM(fct_str(jj))//","//&
3553	& TRIM(fct_str(jk))//")" )
3554	!WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
3555	!WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
3556	ENDIF
3557	! and check it has never gone negative
3558	IF( tg_gdepw_0%d_value(ji,jj,jk,1) < 0._dp .OR. &
3559	& tg_gdept_0%d_value(ji,jj,jk,1) < 0._dp ) THEN
3560	CALL logger_fatal('GRId ZGR SCO : gdepw or gdept =< 0 at point ('//&
3561	& TRIM(fct_str(ji))//","//&
3562	& TRIM(fct_str(jj))//","//&
3563	& TRIM(fct_str(jk))//")" )
3564	!WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
3565	!WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
3566	ENDIF
3567	! and check it never exceeds the total depth
3568	IF( tg_gdepw_0%d_value(ji,jj,jk,1) > tg_hbatt%d_value(ji,jj,1,1) ) THEN
3569	CALL logger_fatal('GRID ZGR SCO: gdepw > hbatt at point ('//&
3570	& TRIM(fct_str(ji))//","//&
3571	& TRIM(fct_str(jj))//","//&
3572	& TRIM(fct_str(jk))//")" )
3573	!WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
3574	ENDIF
3575	ENDDO
3576
3577	DO jk = 1, INT(tg_mbathy%d_value(ji,jj,1,1),i4)-1
3578	! and check it never exceeds the total depth
3579	IF( tg_gdept_0%d_value(ji,jj,jk,1) > tg_hbatt%d_value(ji,jj,1,1) ) THEN
3580	CALL logger_fatal('GRID ZGR SCO: gdept > hbatt at point ('//&
3581	& TRIM(fct_str(ji))//","//&
3582	& TRIM(fct_str(jj))//","//&
3583	& TRIM(fct_str(jk))//")" )
3584	!WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
3585	ENDIF
3586	ENDDO
3587	ENDIF
3588
3589	ENDDO
3590	ENDDO
3591
3592	END SUBROUTINE grid_zgr__sco_fill
3593	!-------------------------------------------------------------------
3594	!> @brief This subroutine stretch the s-coordinate system
3595	!>
3596	!> @details
3597	!> ** Method : s-coordinate stretch using the Song and Haidvogel 1994
3598	!> mixed S/sigma coordinate
3599	!>
3600	!> Reference : Song and Haidvogel 1994.
3601	!>
3602	!> @author J.Paul
3603	!> @date September, 2015 - rewrite from domzgr
3604	!> @date October, 2016
3605	!> - add wetting and drying option
3606	!>
3607	!> @param[in] td_nam
3608	!> @param[in] jpi
3609	!> @param[in] jpj
3610	!> @param[in] jpk
3611	!> @param[in] dd_scosrf
3612	!-------------------------------------------------------------------
3613	SUBROUTINE grid_zgr__sco_s_sh94( td_nam,jpi,jpj,jpk, &
3614	& dd_scosrf )
3615	IMPLICIT NONE
3616	! Argument
3617	TYPE(TNAMZ), INTENT(IN ) :: td_nam
3618	INTEGER(i4), INTENT(IN ) :: jpi
3619	INTEGER(i4), INTENT(IN ) :: jpj
3620	INTEGER(i4), INTENT(IN ) :: jpk
3621	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_scosrf
3622
3623	! local variable
3624	REAL(dp) :: zcoeft
3625	REAL(dp) :: zcoefw
3626
3627	REAL(dp) :: ztmpu
3628	REAL(dp) :: ztmpv
3629	REAL(dp) :: ztmpf
3630	REAL(dp) :: ztmpu1
3631	REAL(dp) :: ztmpv1
3632	REAL(dp) :: ztmpf1
3633
3634	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigw3
3635	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigt3
3636	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsi3w3
3637	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigt3
3638	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigw3
3639	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtu3
3640	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtv3
3641	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtf3
3642	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigwu3
3643	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigwv3
3644
3645	! loop indices
3646	INTEGER(i4) :: ji
3647	INTEGER(i4) :: jj
3648	INTEGER(i4) :: jk
3649	!----------------------------------------------------------------
3650
3651	ALLOCATE( z_gsigw3(jpi,jpj,jpk))
3652	ALLOCATE( z_gsigt3(jpi,jpj,jpk))
3653	ALLOCATE( z_gsi3w3(jpi,jpj,jpk))
3654	ALLOCATE( z_esigt3(jpi,jpj,jpk))
3655	ALLOCATE( z_esigw3(jpi,jpj,jpk))
3656	ALLOCATE( z_esigtu3(jpi,jpj,jpk))
3657	ALLOCATE( z_esigtv3(jpi,jpj,jpk))
3658	ALLOCATE( z_esigtf3(jpi,jpj,jpk))
3659	ALLOCATE( z_esigwu3(jpi,jpj,jpk))
3660	ALLOCATE( z_esigwv3(jpi,jpj,jpk))
3661
3662	z_gsigw3(:,:,:) =0._dp
3663	z_gsigt3(:,:,:) =0._dp
3664	z_gsi3w3(:,:,:) =0._dp
3665	z_esigt3(:,:,:) =0._dp
3666	z_esigw3(:,:,:) =0._dp
3667
3668	z_esigtu3(:,:,:)=0._dp
3669	z_esigtv3(:,:,:)=0._dp
3670	z_esigtf3(:,:,:)=0._dp
3671	z_esigwu3(:,:,:)=0._dp
3672	z_esigwv3(:,:,:)=0._dp
3673
3674	DO ji = 1, jpi
3675	DO jj = 1, jpj
3676
3677	IF( tg_hbatt%d_value(ji,jj,1,1) > td_nam%d_hc ) THEN !deep water, stretched sigma
3678	DO jk = 1, jpk
3679	z_gsigw3(ji,jj,jk) = -grid_zgr__sco_fssig1( td_nam, jpk, REAL(jk,dp)-0.5_dp, td_nam%d_bb )
3680	z_gsigt3(ji,jj,jk) = -grid_zgr__sco_fssig1( td_nam, jpk, REAL(jk,dp) , td_nam%d_bb )
3681	END DO
3682	ELSE ! shallow water, uniform sigma
3683	DO jk = 1, jpk
3684	z_gsigw3(ji,jj,jk) = REAL(jk-1,dp) / REAL(jpk-1,dp)
3685	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,dp) + 0.5_dp ) / REAL(jpk-1,dp)
3686	END DO
3687	ENDIF
3688
3689	DO jk = 1, jpk-1
3690	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
3691	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
3692	END DO
3693	z_esigw3(ji,jj,1 ) = 2._dp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
3694	z_esigt3(ji,jj,jpk) = 2._dp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
3695
3696	! Coefficients for vertical depth as the sum of e3w scale factors
3697	z_gsi3w3(ji,jj,1) = 0.5_dp * z_esigw3(ji,jj,1)
3698	DO jk = 2, jpk
3699	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
3700	END DO
3701	!
3702	DO jk = 1, jpk
3703	zcoeft = ( REAL(jk,dp) - 0.5_dp ) / REAL(jpk-1,dp)
3704	zcoefw = ( REAL(jk,dp) - 1.0_dp ) / REAL(jpk-1,dp)
3705	tg_gdept_0%d_value(ji,jj,jk,1) = ( dd_scosrf(ji,jj) &
3706	& + (tg_hbatt%d_value(ji,jj,1,1)-td_nam%d_hc)*z_gsigt3(ji,jj,jk) &
3707	& + td_nam%d_hc*zcoeft )
3708	tg_gdepw_0%d_value(ji,jj,jk,1) = ( dd_scosrf(ji,jj) &
3709	& + (tg_hbatt%d_value(ji,jj,1,1)-td_nam%d_hc)*z_gsigw3(ji,jj,jk) &
3710	& + td_nam%d_hc*zcoefw )
3711	!tg_gdep3w_0%d_value(ji,jj,jk,1) = ( dd_scosrf(ji,jj) &
3712	! & + (tg_hbatt%d_value(ji,jj,1,1)-td_nam%d_hc)*z_gsi3w3(ji,jj,jk) &
3713	! & + td_nam%d_hc*zcoeft )
3714	END DO
3715
3716	END DO ! for all jj's
3717	END DO ! for all ji's
3718
3719	DO ji = 1, jpi-1
3720	DO jj = 1, jpj-1
3721	! extended for Wetting/Drying case
3722	ztmpu = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji+1,jj ,1,1)
3723	ztmpv = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji ,jj+1,1,1)
3724	ztmpf = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji+1,jj ,1,1) + &
3725	& tg_hbatt%d_value(ji ,jj+1,1,1) + tg_hbatt%d_value(ji+1,jj+1,1,1)
3726
3727	ztmpu1 = tg_hbatt%d_value(ji ,jj ,1,1) * tg_hbatt%d_value(ji+1,jj ,1,1)
3728	ztmpv1 = tg_hbatt%d_value(ji ,jj ,1,1) * tg_hbatt%d_value(ji ,jj+1,1,1)
3729	ztmpf1 = MIN(tg_hbatt%d_value(ji ,jj ,1,1), tg_hbatt%d_value(ji+1,jj ,1,1), &
3730	& tg_hbatt%d_value(ji ,jj+1,1,1), tg_hbatt%d_value(ji+1,jj+1,1,1)) &
3731	& * MAX(tg_hbatt%d_value(ji ,jj ,1,1), tg_hbatt%d_value(ji+1,jj ,1,1), &
3732	& tg_hbatt%d_value(ji ,jj+1,1,1), tg_hbatt%d_value(ji+1,jj+1,1,1))
3733
3734	DO jk = 1, jpk
3735
3736	IF( td_nam%l_wd .AND. &
3737	& ( ztmpu1 < 0._dp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
3738	z_esigtu3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
3739	z_esigwu3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
3740	ELSE
3741	z_esigtu3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj,1,1)*z_esigt3(ji ,jj,jk) &
3742	& + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigt3(ji+1,jj,jk) ) / ztmpu
3743	z_esigwu3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj,1,1)*z_esigw3(ji ,jj,jk) &
3744	& + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigw3(ji+1,jj,jk)) / ztmpu
3745	ENDIF
3746
3747	IF( td_nam%l_wd .AND. &
3748	& ( ztmpv1 < 0._dp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
3749	z_esigtv3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
3750	z_esigwv3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
3751	ELSE
3752	z_esigtv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj ,1,1)*z_esigt3(ji,jj ,jk) &
3753	& + tg_hbatt%d_value(ji,jj+1,1,1)*z_esigt3(ji,jj+1,jk)) / ztmpv
3754	z_esigwv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj ,1,1)*z_esigw3(ji,jj ,jk) &
3755	& + tg_hbatt%d_value(ji,jj+1,1,1)*z_esigw3(ji,jj+1,jk)) / ztmpv
3756	ENDIF
3757
3758	IF( td_nam%l_wd .AND. &
3759	& ( ztmpf1 < 0._dp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
3760	z_esigtf3(ji,jj,jk) = 0.25_dp * ( z_esigt3(ji ,jj ,jk) &
3761	& + z_esigt3(ji+1,jj ,jk) &
3762	& + z_esigt3(ji ,jj+1,jk) &
3763	& + z_esigt3(ji+1,jj+1,jk) )
3764	ELSE
3765	z_esigtf3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj ,1,1)*z_esigt3(ji ,jj ,jk) &
3766	& + tg_hbatt%d_value(ji+1,jj ,1,1)*z_esigt3(ji+1,jj ,jk) &
3767	& + tg_hbatt%d_value(ji ,jj+1,1,1)*z_esigt3(ji ,jj+1,jk) &
3768	& + tg_hbatt%d_value(ji+1,jj+1,1,1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
3769	ENDIF
3770	!
3771	tg_e3t_0%d_value(ji,jj,jk,1) = ( ( tg_hbatt%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigt3 (ji,jj,jk) &
3772	& + td_nam%d_hc / REAL(jpk-1,dp) )
3773	tg_e3u_0%d_value(ji,jj,jk,1) = ( ( tg_hbatu%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigtu3(ji,jj,jk) &
3774	& + td_nam%d_hc / REAL(jpk-1,dp) )
3775	tg_e3v_0%d_value(ji,jj,jk,1) = ( ( tg_hbatv%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigtv3(ji,jj,jk) &
3776	& + td_nam%d_hc / REAL(jpk-1,dp) )
3777	tg_e3f_0%d_value(ji,jj,jk,1) = ( ( tg_hbatf%d_value(ji,jj,1,1) - td_nam%d_hc ) *z_esigtf3(ji,jj,jk) &
3778	& + td_nam%d_hc/REAL(jpk-1,dp) )
3779
3780	tg_e3w_0%d_value (ji,jj,jk,1)= ( ( tg_hbatt%d_value(ji,jj,1,1) - td_nam%d_hc )*z_esigw3 (ji,jj,jk) &
3781	& + td_nam%d_hc / REAL(jpk-1,dp) )
3782	tg_e3uw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatu%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwu3(ji,jj,jk) &
3783	& + td_nam%d_hc/REAL(jpk-1,dp) )
3784	tg_e3vw_0%d_value(ji,jj,jk,1)= ( ( tg_hbatv%d_value(ji,jj,1,1) - td_nam%d_hc)*z_esigwv3(ji,jj,jk) &
3785	& + td_nam%d_hc/REAL(jpk-1,dp) )
3786	END DO
3787	END DO
3788	END DO
3789
3790	DEALLOCATE( z_gsigw3 )
3791	DEALLOCATE( z_gsigt3 )
3792	DEALLOCATE( z_gsi3w3 )
3793	DEALLOCATE( z_esigt3 )
3794	DEALLOCATE( z_esigw3 )
3795	DEALLOCATE( z_esigtu3 )
3796	DEALLOCATE( z_esigtv3 )
3797	DEALLOCATE( z_esigtf3 )
3798	DEALLOCATE( z_esigwu3 )
3799	DEALLOCATE( z_esigwv3 )
3800
3801	END SUBROUTINE grid_zgr__sco_s_sh94
3802	!-------------------------------------------------------------------
3803	!> @brief This subroutine stretch the s-coordinate system
3804	!>
3805	!> ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
3806	!> mixed S/sigma/Z coordinate
3807	!>
3808	!> This method allows the maintenance of fixed surface and or
3809	!> bottom cell resolutions (cf. geopotential coordinates)
3810	!> within an analytically derived stretched S-coordinate framework.
3811	!>
3812	!>
3813	!> Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
3814	!>
3815	!> @author J.Paul
3816	!> @date September, 2015 - rewrite from domzgr
3817	!> @date October, 2016
3818	!> - add wetting and drying option
3819	!>
3820	!> @param[in] td_nam
3821	!> @param[in] jpi
3822	!> @param[in] jpj
3823	!> @param[in] jpk
3824	!> @param[in] dd_scosrf
3825	!-------------------------------------------------------------------
3826	SUBROUTINE grid_zgr__sco_s_sf12( td_nam,jpi,jpj,jpk, &
3827	& dd_scosrf )
3828	IMPLICIT NONE
3829	! Argument
3830	TYPE(TNAMZ), INTENT(IN ) :: td_nam
3831	INTEGER(i4), INTENT(IN ) :: jpi
3832	INTEGER(i4), INTENT(IN ) :: jpj
3833	INTEGER(i4), INTENT(IN ) :: jpk
3834	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_scosrf
3835
3836	! local variable
3837	REAL(dp) :: zsmth ! smoothing around critical depth
3838	REAL(dp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
3839
3840	REAL(dp) :: ztmpu
3841	REAL(dp) :: ztmpv
3842	REAL(dp) :: ztmpf
3843	REAL(dp) :: ztmpu1
3844	REAL(dp) :: ztmpv1
3845	REAL(dp) :: ztmpf1
3846
3847	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigw3
3848	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsigt3
3849	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_gsi3w3
3850	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigt3
3851	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigw3
3852	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtu3
3853	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtv3
3854	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigtf3
3855	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigwu3
3856	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: z_esigwv3
3857
3858	! loop indices
3859	INTEGER(i4) :: ji
3860	INTEGER(i4) :: jj
3861	INTEGER(i4) :: jk
3862	!----------------------------------------------------------------
3863	ALLOCATE( z_gsigw3(jpi,jpj,jpk))
3864	ALLOCATE( z_gsigt3(jpi,jpj,jpk))
3865	ALLOCATE( z_gsi3w3(jpi,jpj,jpk))
3866	ALLOCATE( z_esigt3(jpi,jpj,jpk))
3867	ALLOCATE( z_esigw3(jpi,jpj,jpk))
3868	ALLOCATE( z_esigtu3(jpi,jpj,jpk))
3869	ALLOCATE( z_esigtv3(jpi,jpj,jpk))
3870	ALLOCATE( z_esigtf3(jpi,jpj,jpk))
3871	ALLOCATE( z_esigwu3(jpi,jpj,jpk))
3872	ALLOCATE( z_esigwv3(jpi,jpj,jpk))
3873
3874	z_gsigw3(:,:,:) =0._dp
3875	z_gsigt3(:,:,:) =0._dp
3876	z_gsi3w3(:,:,:) =0._dp
3877	z_esigt3(:,:,:) =0._dp
3878	z_esigw3(:,:,:) =0._dp
3879	z_esigtu3(:,:,:)=0._dp
3880	z_esigtv3(:,:,:)=0._dp
3881	z_esigtf3(:,:,:)=0._dp
3882	z_esigwu3(:,:,:)=0._dp
3883	z_esigwv3(:,:,:)=0._dp
3884
3885	DO ji = 1, jpi
3886	DO jj = 1, jpj
3887
3888	IF( tg_hbatt%d_value(ji,jj,1,1) > td_nam%d_hc )THEN !deep water, stretched sigma
3889
3890	! this forces a linear bottom cell depth relationship with H,.
3891	! could be changed by users but care must be taken to do so carefully
3892	zzb = tg_hbatt%d_value(ji,jj,1,1)*td_nam%d_zb_a + td_nam%d_zb_b
3893
3894	zzb = 1.0_dp-(zzb/tg_hbatt%d_value(ji,jj,1,1))
3895
3896	zzs = td_nam%d_zs / tg_hbatt%d_value(ji,jj,1,1)
3897
3898	IF( td_nam%d_efold /= 0.0_dp )THEN
3899	zsmth = TANH( (tg_hbatt%d_value(ji,jj,1,1)-td_nam%d_hc ) / td_nam%d_efold )
3900	ELSE
3901	zsmth = 1.0_dp
3902	ENDIF
3903
3904	DO jk = 1, jpk
3905	z_gsigw3(ji,jj,jk) = REAL(jk-1,dp) /REAL(jpk-1,dp)
3906	z_gsigt3(ji,jj,jk) = (REAL(jk-1,dp)+0.5_dp)/REAL(jpk-1,dp)
3907	ENDDO
3908	z_gsigw3(ji,jj,:) = grid_zgr__sco_fgamma( td_nam, jpk, z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
3909	z_gsigt3(ji,jj,:) = grid_zgr__sco_fgamma( td_nam, jpk, z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
3910
3911	ELSE IF( td_nam%l_sigcrit )THEN ! shallow water, uniform sigma
3912
3913	DO jk = 1, jpk
3914	z_gsigw3(ji,jj,jk) = REAL(jk-1,dp) /REAL(jpk-1,dp)
3915	z_gsigt3(ji,jj,jk) = (REAL(jk-1,dp)+0.5)/REAL(jpk-1,dp)
3916	END DO
3917
3918	ELSE ! shallow water, z coordinates
3919
3920	DO jk = 1, jpk
3921	z_gsigw3(ji,jj,jk) = REAL(jk-1,dp) /REAL(jpk-1,dp)*(td_nam%d_hc/tg_hbatt%d_value(ji,jj,1,1))
3922	z_gsigt3(ji,jj,jk) = (REAL(jk-1,dp)+0.5_dp)/REAL(jpk-1,dp)*(td_nam%d_hc/tg_hbatt%d_value(ji,jj,1,1))
3923	END DO
3924
3925	ENDIF
3926
3927	DO jk = 1, jpk-1
3928	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
3929	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
3930	END DO
3931	z_esigw3(ji,jj,1 ) = 2.0_dp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
3932	z_esigt3(ji,jj,jpk) = 2.0_dp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
3933
3934	! Coefficients for vertical depth as the sum of e3w scale factors
3935	z_gsi3w3(ji,jj,1) = 0.5_dp * z_esigw3(ji,jj,1)
3936	DO jk = 2, jpk
3937	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
3938	END DO
3939
3940	DO jk = 1, jpk
3941	tg_gdept_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigt3(ji,jj,jk)
3942	tg_gdepw_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsigw3(ji,jj,jk)
3943	!tg_gdep3w_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_gsi3w3(ji,jj,jk)
3944	END DO
3945
3946	ENDDO ! for all jj's
3947	ENDDO ! for all ji's
3948
3949	DO ji=1,jpi-1
3950	DO jj=1,jpj-1
3951
3952	! extend to suit for Wetting/Drying case
3953	ztmpu = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji+1,jj ,1,1)
3954	ztmpv = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji ,jj+1,1,1)
3955	ztmpf = tg_hbatt%d_value(ji ,jj ,1,1) + tg_hbatt%d_value(ji+1,jj ,1,1) + &
3956	& tg_hbatt%d_value(ji ,jj+1,1,1) + tg_hbatt%d_value(ji+1,jj+1,1,1)
3957
3958	ztmpu1 = tg_hbatt%d_value(ji ,jj ,1,1) * tg_hbatt%d_value(ji+1,jj ,1,1)
3959	ztmpv1 = tg_hbatt%d_value(ji ,jj ,1,1) * tg_hbatt%d_value(ji ,jj+1,1,1)
3960	ztmpf1 = MIN(tg_hbatt%d_value(ji ,jj ,1,1), tg_hbatt%d_value(ji+1,jj ,1,1), &
3961	& tg_hbatt%d_value(ji ,jj+1,1,1), tg_hbatt%d_value(ji+1,jj+1,1,1)) &
3962	& * MAX(tg_hbatt%d_value(ji ,jj ,1,1), tg_hbatt%d_value(ji+1,jj ,1,1), &
3963	& tg_hbatt%d_value(ji ,jj+1,1,1), tg_hbatt%d_value(ji+1,jj+1,1,1))
3964
3965	DO jk = 1, jpk
3966
3967	IF( td_nam%l_wd .AND. &
3968	& ( ztmpu1 < 0._dp .OR. ABS(ztmpu) < td_nam%d_wdmin1 ) )THEN
3969	z_esigtu3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
3970	z_esigwu3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
3971	ELSE
3972	z_esigtu3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj,1,1)*z_esigt3(ji ,jj,jk) &
3973	& + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigt3(ji+1,jj,jk) ) / ztmpu
3974	z_esigwu3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj,1,1)*z_esigw3(ji ,jj,jk) &
3975	& + tg_hbatt%d_value(ji+1,jj,1,1)*z_esigw3(ji+1,jj,jk) ) / ztmpu
3976	ENDIF
3977
3978	IF( td_nam%l_wd .AND. &
3979	& ( ztmpv1 < 0._dp .OR. ABS(ztmpv) < td_nam%d_wdmin1 ) )THEN
3980	z_esigtv3(ji,jj,jk) = 0.5_dp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
3981	z_esigwv3(ji,jj,jk) = 0.5_dp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
3982	ELSE
3983	z_esigtv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj ,1,1)*z_esigt3(ji,jj ,jk) &
3984	& + tg_hbatt%d_value(ji,jj+1,1,1)*z_esigt3(ji,jj+1,jk)) / ztmpv
3985	z_esigwv3(ji,jj,jk) = ( tg_hbatt%d_value(ji,jj ,1,1)*z_esigw3(ji,jj ,jk) &
3986	& + tg_hbatt%d_value(ji,jj+1,1,1)*z_esigw3(ji,jj+1,jk)) / ztmpv
3987	ENDIF
3988
3989	IF( td_nam%l_wd .AND. &
3990	& ( ztmpf1 < 0._dp .OR. ABS(ztmpf) < td_nam%d_wdmin1 ) )THEN
3991	z_esigtf3(ji,jj,jk) = 0.25_dp * ( z_esigt3(ji ,jj ,jk) &
3992	& + z_esigt3(ji+1,jj ,jk) &
3993	& + z_esigt3(ji ,jj+1,jk) &
3994	& + z_esigt3(ji+1,jj+1,jk) )
3995	ELSE
3996	z_esigtf3(ji,jj,jk) = ( tg_hbatt%d_value(ji ,jj ,1,1)*z_esigt3(ji ,jj ,jk) &
3997	& + tg_hbatt%d_value(ji+1,jj ,1,1)*z_esigt3(ji+1,jj ,jk) &
3998	& + tg_hbatt%d_value(ji ,jj+1,1,1)*z_esigt3(ji ,jj+1,jk) &
3999	& + tg_hbatt%d_value(ji+1,jj+1,1,1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
4000	ENDIF
4001
4002	tg_e3t_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatt%d_value(ji,jj,1,1))*z_esigt3 (ji,jj,jk)
4003	tg_e3u_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatu%d_value(ji,jj,1,1))*z_esigtu3(ji,jj,jk)
4004	tg_e3v_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatv%d_value(ji,jj,1,1))*z_esigtv3(ji,jj,jk)
4005	tg_e3f_0%d_value(ji,jj,jk,1) = (dd_scosrf(ji,jj)+tg_hbatf%d_value(ji,jj,1,1))*z_esigtf3(ji,jj,jk)
4006	!
4007	tg_e3w_0%d_value(ji,jj,jk,1) = tg_hbatt%d_value(ji,jj,1,1)*z_esigw3 (ji,jj,jk)
4008	tg_e3uw_0%d_value(ji,jj,jk,1) = tg_hbatu%d_value(ji,jj,1,1)*z_esigwu3(ji,jj,jk)
4009	tg_e3vw_0%d_value(ji,jj,jk,1) = tg_hbatv%d_value(ji,jj,1,1)*z_esigwv3(ji,jj,jk)
4010	END DO
4011
4012	ENDDO
4013	ENDDO
4014
4015	CALL lbc_lnk(tg_e3t_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4016	CALL lbc_lnk(tg_e3u_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4017	CALL lbc_lnk(tg_e3v_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4018	CALL lbc_lnk(tg_e3f_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4019	CALL lbc_lnk(tg_e3w_0 %d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4020	CALL lbc_lnk(tg_e3uw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4021	CALL lbc_lnk(tg_e3vw_0%d_value(:,:,:,1),'T', td_nam%i_perio, 1._dp)
4022
4023	DEALLOCATE( z_gsigw3 )
4024	DEALLOCATE( z_gsigt3 )
4025	DEALLOCATE( z_gsi3w3 )
4026	DEALLOCATE( z_esigt3 )
4027	DEALLOCATE( z_esigw3 )
4028	DEALLOCATE( z_esigtu3 )
4029	DEALLOCATE( z_esigtv3 )
4030	DEALLOCATE( z_esigtf3 )
4031	DEALLOCATE( z_esigwu3 )
4032	DEALLOCATE( z_esigwv3 )
4033
4034	END SUBROUTINE grid_zgr__sco_s_sf12
4035	!-------------------------------------------------------------------
4036	!> @brief This subroutine stretch the s-coordinate system
4037	!>
4038	!>
4039	!> ** Method : s-coordinate stretch
4040	!>
4041	!> Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
4042	!>
4043	!> @author J.Paul
4044	!> @date September, 2015 - rewrite from domzgr
4045	!>
4046	!> @param[in] td_nam
4047	!> @param[in] jpi
4048	!> @param[in] jpj
4049	!> @param[in] jpk
4050	!> @param[in] dd_scosrf
4051	!> @param[in] dd_hift
4052	!> @param[in] dd_hifu
4053	!> @param[in] dd_hifv
4054	! @param[in] dd_hiff
4055	!-------------------------------------------------------------------
4056	SUBROUTINE grid_zgr__sco_s_tanh( td_nam,jpi,jpj,jpk, &
4057	& dd_scosrf, &
4058	& dd_hift, dd_hifu, dd_hifv, dd_hiff )
4059	IMPLICIT NONE
4060	! Argument
4061	TYPE(TNAMZ), INTENT(IN ) :: td_nam
4062	INTEGER(i4), INTENT(IN ) :: jpi
4063	INTEGER(i4), INTENT(IN ) :: jpj
4064	INTEGER(i4), INTENT(IN ) :: jpk
4065	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_scosrf
4066	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_hift
4067	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_hifu
4068	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_hifv
4069	REAL(dp) , DIMENSION(:,:), INTENT(IN ) :: dd_hiff
4070
4071	! local variable
4072	REAL(dp) :: zcoeft
4073	REAL(dp) :: zcoefw
4074
4075	! loop indices
4076	INTEGER(i4) :: ji
4077	INTEGER(i4) :: jj
4078	INTEGER(i4) :: jk
4079	!----------------------------------------------------------------
4080
4081	tg_gsigt%d_value(1,1,:,1) =0._dp
4082	tg_gsigw%d_value(1,1,:,1) =0._dp
4083	tg_gsi3w%d_value(1,1,:,1) =0._dp
4084	tg_esigt%d_value(1,1,:,1) =0._dp
4085	tg_esigw%d_value(1,1,:,1) =0._dp
4086
4087	DO jk = 1, jpk
4088	tg_gsigw%d_value(1,1,jk,1) = -grid_zgr__sco_fssig( td_nam, jpk, REAL(jk,dp)-0.5_dp )
4089	tg_gsigt%d_value(1,1,jk,1) = -grid_zgr__sco_fssig( td_nam, jpk, REAL(jk,dp) )
4090	END DO
4091	CALL logger_info('z_gsigw 1 jpk '//TRIM(fct_str(tg_gsigw%d_value(1,1,1,1)))//&
4092	& TRIM(fct_str(tg_gsigw%d_value(1,1,jpk,1))) )
4093
4094	! Coefficients for vertical scale factors at w-, t- levels
4095	!!gm bug : define it from analytical function, not like juste bellow....
4096	!!gm or betteroffer the 2 possibilities....
4097	DO jk = 1, jpk-1
4098	tg_esigt%d_value(1,1,jk ,1) = tg_gsigw%d_value(1,1,jk+1,1) - tg_gsigw%d_value(1,1,jk,1)
4099	tg_esigw%d_value(1,1,jk+1,1) = tg_gsigt%d_value(1,1,jk+1,1) - tg_gsigt%d_value(1,1,jk,1)
4100	END DO
4101	tg_esigw%d_value(1,1, 1 ,1) = 2._dp * ( tg_gsigt%d_value(1,1,1 ,1) - tg_gsigw%d_value(1,1,1 ,1) )
4102	tg_esigt%d_value(1,1,jpk,1) = 2._dp * ( tg_gsigt%d_value(1,1,jpk,1) - tg_gsigw%d_value(1,1,jpk,1) )
4103
4104	! Coefficients for vertical depth as the sum of e3w scale factors
4105	tg_gsi3w%d_value(1,1,1,1) = 0.5_dp * tg_esigw%d_value(1,1,1,1)
4106	DO jk = 2, jpk
4107	tg_gsi3w%d_value(1,1,jk,1) = tg_gsi3w%d_value(1,1,jk-1,1) + tg_esigw%d_value(1,1,jk,1)
4108	END DO
4109	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
4110	DO jk = 1, jpk
4111	zcoeft = ( REAL(jk,dp) - 0.5_dp ) / REAL(jpk-1,dp)
4112	zcoefw = ( REAL(jk,dp) - 1.0_dp ) / REAL(jpk-1,dp)
4113	tg_gdept_0%d_value (:,:,jk,1) = dd_scosrf(:,:) + ( tg_hbatt%d_value(:,:,1 ,1) - dd_hift(:,:) ) &
4114	& * tg_gsigt%d_value(1,1,jk,1) + dd_hift(:,:)*zcoeft
4115	tg_gdepw_0%d_value (:,:,jk,1) = dd_scosrf(:,:) + ( tg_hbatt%d_value(:,:,1 ,1) - dd_hift(:,:) ) &
4116	& * tg_gsigw%d_value(1,1,jk,1) + dd_hift(:,:)*zcoefw
4117	!tg_gdep3w_0%d_value(:,:,jk,1) = dd_scosrf(:,:) + ( tg_hbatt%d_value(:,:,1 ,1) - dd_hift(:,:)) &
4118	! * tg_gsi3w%d_value(1,1,jk,1) + dd_hift(:,:)*zcoeft
4119	END DO
4120	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
4121	DO jj = 1, jpj
4122	DO ji = 1, jpi
4123	DO jk = 1, jpk
4124	tg_e3t_0%d_value(ji,jj,jk,1) = ( (tg_hbatt%d_value(ji,jj,1 ,1) - dd_hift(ji,jj)) &
4125	& * tg_esigt%d_value(1 ,1 ,jk,1) + dd_hift(ji,jj)/REAL(jpk-1,dp) )
4126	tg_e3u_0%d_value(ji,jj,jk,1) = ( (tg_hbatu%d_value(ji,jj,1 ,1) - dd_hifu(ji,jj)) &
4127	& * tg_esigt%d_value(1 ,1 ,jk,1) + dd_hifu(ji,jj)/REAL(jpk-1,dp) )
4128	tg_e3v_0%d_value(ji,jj,jk,1) = ( (tg_hbatv%d_value(ji,jj,1 ,1) - dd_hifv(ji,jj)) &
4129	& * tg_esigt%d_value(1 ,1 ,jk,1) + dd_hifv(ji,jj)/REAL(jpk-1,dp) )
4130	tg_e3f_0%d_value(ji,jj,jk,1) = ( (tg_hbatf%d_value(ji,jj,1 ,1) - dd_hiff(ji,jj)) &
4131	& * tg_esigt%d_value(1 ,1, jk,1) + dd_hiff(ji,jj)/REAL(jpk-1,dp) )
4132
4133	tg_e3w_0%d_value (ji,jj,jk,1)= ( (tg_hbatt%d_value(ji,jj,1 ,1) - dd_hift(ji,jj)) &
4134	& * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hift(ji,jj)/REAL(jpk-1,dp) )
4135	tg_e3uw_0%d_value(ji,jj,jk,1)= ( (tg_hbatu%d_value(ji,jj,1 ,1) - dd_hifu(ji,jj)) &
4136	& * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifu(ji,jj)/REAL(jpk-1,dp) )
4137	tg_e3vw_0%d_value(ji,jj,jk,1)= ( (tg_hbatv%d_value(ji,jj,1 ,1) - dd_hifv(ji,jj)) &
4138	& * tg_esigw%d_value(1 ,1 ,jk,1) + dd_hifv(ji,jj)/REAL(jpk-1,dp) )
4139	END DO
4140	END DO
4141	END DO
4142
4143	END SUBROUTINE grid_zgr__sco_s_tanh
4144	!!----------------------------------------------------------------------
4145	!> @brief This function provide the analytical function in s-coordinate
4146	!>
4147	!> @details
4148	!> ** Method : the function provide the non-dimensional position of
4149	!> T and W (i.e. between 0 and 1)
4150	!> T-points at integer values (between 1 and jpk)
4151	!> W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
4152	!>
4153	!> @author J.Paul
4154	!> @date September, 2015 - rewrite from domzgr
4155	!>
4156	!> @param[in] td_nam
4157	!> @param[in] jpk
4158	!> @param[in] pk
4159	!!----------------------------------------------------------------------
4160	FUNCTION grid_zgr__sco_fssig( td_nam, jpk, pk ) RESULT( pf )
4161	IMPLICIT NONE
4162	! Argument
4163	TYPE(TNAMZ), INTENT(IN ) :: td_nam
4164	INTEGER(i4), INTENT(IN ) :: jpk
4165	REAL(dp) , INTENT(IN ) :: pk ! continuous "k" coordinate
4166	REAL(dp) :: pf ! sigma value
4167	!!----------------------------------------------------------------------
4168	!
4169	pf = ( TANH( td_nam%d_theta * ( -(pk-0.5_dp) / REAL(jpk-1,dp) + td_nam%d_thetb ) ) &
4170	& - TANH( td_nam%d_thetb * td_nam%d_theta ) ) &
4171	& * ( COSH( td_nam%d_theta ) &
4172	& + COSH( td_nam%d_theta * ( 2._dp * td_nam%d_thetb - 1._dp ) ) ) &
4173	& / ( 2._dp * SINH( td_nam%d_theta ) )
4174	!
4175	END FUNCTION grid_zgr__sco_fssig
4176	!!----------------------------------------------------------------------
4177	!> @brief This function provide the Song and Haidvogel version of the analytical function in s-coordinate
4178	!>
4179	!> @details
4180	!> ** Method : the function provides the non-dimensional position of
4181	!> T and W (i.e. between 0 and 1)
4182	!> T-points at integer values (between 1 and jpk)
4183	!> W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
4184	!>
4185	!> @author J.Paul
4186	!> @date September, 2015 - rewrite from domzgr
4187	!>
4188	!> @param[in] td_nam
4189	!> @param[in] jpi
4190	!> @param[in] jpj
4191	!> @param[in] pk1
4192	!> @param[in] pbb
4193	!!----------------------------------------------------------------------
4194	FUNCTION grid_zgr__sco_fssig1( td_nam, jpk, pk1, pbb ) RESULT( pf1 )
4195	IMPLICIT NONE
4196	! Argument
4197	TYPE(TNAMZ), INTENT(IN ) :: td_nam
4198	INTEGER(i4), INTENT(IN ) :: jpk
4199	REAL(dp) , INTENT(IN ) :: pk1 ! continuous "k" coordinate
4200	REAL(dp) , INTENT(IN ) :: pbb ! Stretching coefficient
4201	REAL(dp) :: pf1 ! sigma value
4202	!!----------------------------------------------------------------------
4203	!
4204	IF ( td_nam%d_theta == 0 ) then ! uniform sigma
4205	pf1 = - ( pk1 - 0.5_dp ) / REAL(jpk-1,dp)
4206	ELSE ! stretched sigma
4207	pf1 = ( 1._dp - pbb ) * ( SINH( td_nam%d_theta*(-(pk1-0.5_dp)/REAL(jpk-1,dp)) ) ) / SINH( td_nam%d_theta ) &
4208	& + pbb * ( (TANH( td_nam%d_theta( (-(pk1-0.5_dp)/REAL(jpk-1,dp)) + 0.5_dp) ) - TANH( 0.5_dp td_nam%d_theta ) ) &
4209	& / ( 2._dp * TANH( 0.5_dp * td_nam%d_theta ) ) )
4210	ENDIF
4211	!
4212	END FUNCTION grid_zgr__sco_fssig1
4213	!!----------------------------------------------------------------------
4214	!> @brief This function provide analytical function for the s-coordinate
4215	!>
4216	!> @details
4217	!> ** Method : the function provides the non-dimensional position of
4218	!> T and W (i.e. between 0 and 1)
4219	!> T-points at integer values (between 1 and jpk)
4220	!> W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
4221	!>
4222	!> This method allows the maintenance of fixed surface and or
4223	!> bottom cell resolutions (cf. geopotential coordinates)
4224	!> within an analytically derived stretched S-coordinate framework.
4225	!>
4226	!> Reference : Siddorn and Furner, in prep
4227	!>
4228	!> @author J.Paul
4229	!> @date September, 2015 - rewrite from domzgr
4230	!>
4231	!> @param[in] td_nam
4232	!> @param[in] jpk
4233	!> @param[in] pk1
4234	!> @param[in] pzb
4235	!> @param[in] pzs
4236	!> @param[in] pzsmth
4237	!!----------------------------------------------------------------------
4238	FUNCTION grid_zgr__sco_fgamma( td_nam, jpk, pk1, pzb, pzs, psmth) RESULT( p_gamma )
4239	IMPLICIT NONE
4240	! Argument
4241	TYPE(TNAMZ), INTENT(IN ) :: td_nam
4242	INTEGER(i4), INTENT(IN ) :: jpk
4243	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: pk1 ! continuous "k" coordinate
4244	! function
4245	REAL(dp) , DIMENSION(jpk) :: p_gamma ! stretched coordinate
4246	! local variable
4247	REAL(dp) , INTENT(IN ) :: pzb ! Bottom box depth
4248	REAL(dp) , INTENT(IN ) :: pzs ! surface box depth
4249	REAL(dp) , INTENT(IN ) :: psmth ! Smoothing parameter
4250	REAL(dp) :: za1,za2,za3 ! local variables
4251	REAL(dp) :: zn1,zn2 ! local variables
4252	REAL(dp) :: za,zb,zx ! local variables
4253
4254	! loop indices
4255	INTEGER(i4) :: jk
4256	!!----------------------------------------------------------------------
4257	!
4258	zn1 = 1._dp / REAL(jpk-1,dp)
4259	zn2 = 1._dp - zn1
4260
4261	za1 = (td_nam%d_alpha+2.0_dp)zn1(td_nam%d_alpha+1.0_dp)-(td_nam%d_alpha+1.0_dp)zn1**(td_nam%d_alpha+2.0_dp)
4262	za2 = (td_nam%d_alpha+2.0_dp)zn2(td_nam%d_alpha+1.0_dp)-(td_nam%d_alpha+1.0_dp)zn2**(td_nam%d_alpha+2.0_dp)
4263	za3 = (zn23.0_dp - za2)/( zn13.0_dp - za1)
4264
4265	za = pzb - za3*(pzs-za1)-za2
4266	za = za/( zn2-0.5_dp(za2+zn22.0_dp) - za3(zn1-0.5_dp(za1+zn1*2.0_dp) ) )
4267	zb = (pzs - za1 - za( zn1-0.5_dp(za1+zn12.0_dp ) ) ) / (zn13.0_dp - za1)
4268	zx = 1.0_dp-za/2.0_dp-zb
4269
4270	DO jk = 1, jpk
4271	p_gamma(jk) = za(pk1(jk)(1.0_dp-pk1(jk)/2.0_dp))+zbpk1(jk)*3.0_dp + &
4272	& zx( (td_nam%d_alpha+2.0_dp)pk1(jk)**(td_nam%d_alpha+1.0_dp)- &
4273	& (td_nam%d_alpha+1.0_dp)pk1(jk)*(td_nam%d_alpha+2.0_dp) )
4274	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_dp-psmth)
4275	ENDDO
4276
4277	!
4278	END FUNCTION grid_zgr__sco_fgamma
4279	!-------------------------------------------------------------------
4280	!> @brief This subroutine stretch the s-coordinate system
4281	!>
4282	!>
4283	!> ** Method : s-coordinate stretch
4284	!>
4285	!> Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
4286	!>
4287	!> @author J.Paul
4288	!> @date September, 2015 - rewrite from domain (dom_stiff)
4289	!>
4290	!> @param[in] td_nam
4291	!> @param[in] jpi
4292	!> @param[in] jpj
4293	!> @param[in] jpk
4294	!-------------------------------------------------------------------
4295	SUBROUTINE grid_zgr_sco_stiff(td_nam, jpi,jpj,jpk)
4296	IMPLICIT NONE
4297	! Argument
4298	TYPE(TNAMZ), INTENT(IN ) :: td_nam
4299	INTEGER(i4), INTENT(IN ) :: jpi
4300	INTEGER(i4), INTENT(IN ) :: jpj
4301	INTEGER(i4), INTENT(IN ) :: jpk
4302
4303	! local variable
4304	REAL(dp) :: zrxmax
4305	REAL(dp), DIMENSION(4) :: zr1
4306
4307	! loop indices
4308	INTEGER(i4) :: ji
4309	INTEGER(i4) :: jj
4310	INTEGER(i4) :: jk
4311	!----------------------------------------------------------------
4312	tg_rx1%d_value(:,:,1,1) = 0._dp
4313
4314	zrxmax = 0._dp
4315	zr1(:) = 0._dp
4316
4317	DO ji = 2, jpi-1
4318	DO jj = 2, jpj-1
4319	DO jk = 1, jpk-1
4320	zr1(1) = tg_umask%d_value(ji-1,jj ,jk,1) &
4321	& * ABS( ( tg_gdepw_0%d_value(ji ,jj ,jk ,1) &
4322	& - tg_gdepw_0%d_value(ji-1,jj ,jk ,1) &
4323	& + tg_gdepw_0%d_value(ji ,jj ,jk+1,1) &
4324	& - tg_gdepw_0%d_value(ji-1,jj ,jk+1,1) ) &
4325	& / ( tg_gdepw_0%d_value(ji ,jj ,jk ,1) &
4326	& + tg_gdepw_0%d_value(ji-1,jj ,jk ,1) &
4327	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1) &
4328	& - tg_gdepw_0%d_value(ji-1,jj ,jk+1,1) &
4329	& + dp_eps) )
4330	zr1(2) = tg_umask%d_value(ji ,jj ,jk,1) &
4331	& * ABS( ( tg_gdepw_0%d_value(ji+1,jj ,jk ,1) &
4332	& - tg_gdepw_0%d_value(ji ,jj ,jk ,1) &
4333	& + tg_gdepw_0%d_value(ji+1,jj ,jk+1,1) &
4334	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1) ) &
4335	& / ( tg_gdepw_0%d_value(ji+1,jj ,jk ,1)&
4336	& + tg_gdepw_0%d_value(ji ,jj ,jk ,1)&
4337	& - tg_gdepw_0%d_value(ji+1,jj ,jk+1,1)&
4338	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1)&
4339	& + dp_eps) )
4340	zr1(3) = tg_vmask%d_value(ji ,jj ,jk,1) &
4341	& * ABS( ( tg_gdepw_0%d_value(ji ,jj+1,jk ,1)&
4342	& - tg_gdepw_0%d_value(ji ,jj ,jk ,1)&
4343	& + tg_gdepw_0%d_value(ji ,jj+1,jk+1,1)&
4344	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1) ) &
4345	& / ( tg_gdepw_0%d_value(ji ,jj+1,jk ,1)&
4346	& + tg_gdepw_0%d_value(ji ,jj ,jk ,1)&
4347	& - tg_gdepw_0%d_value(ji ,jj+1,jk+1,1)&
4348	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1)&
4349	& + dp_eps) )
4350	zr1(4) = tg_vmask%d_value(ji ,jj-1,jk,1) &
4351	& * ABS( ( tg_gdepw_0%d_value(ji ,jj ,jk ,1)&
4352	& - tg_gdepw_0%d_value(ji ,jj-1,jk ,1)&
4353	& + tg_gdepw_0%d_value(ji ,jj ,jk+1,1)&
4354	& - tg_gdepw_0%d_value(ji ,jj-1,jk+1,1) ) &
4355	& / ( tg_gdepw_0%d_value(ji ,jj ,jk ,1)&
4356	& + tg_gdepw_0%d_value(ji ,jj-1,jk ,1)&
4357	& - tg_gdepw_0%d_value(ji ,jj ,jk+1,1)&
4358	& - tg_gdepw_0%d_value(ji ,jj-1,jk+1,1)&
4359	& + dp_eps) )
4360	zrxmax = MAXVAL(zr1(1:4))
4361	tg_rx1%d_value(ji,jj,1,1) = MAX(tg_rx1%d_value(ji,jj,1,1), zrxmax)
4362	END DO
4363	END DO
4364	END DO
4365
4366	CALL lbc_lnk( tg_rx1%d_value(:,:,1,1), 'T', td_nam%i_perio, 1._dp )
4367
4368	zrxmax = MAXVAL(tg_rx1%d_value(:,:,1,1))
4369	CALL logger_info(' GRID ZGR SCO STIFF: maximum grid stiffness ratio: '//&
4370	& TRIM(fct_str(zrxmax)) )
4371
4372	END SUBROUTINE grid_zgr_sco_stiff
4373	END MODULE grid_zgr

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

Download in other formats: