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 branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2016/dev_r6999_CONFIGMAN_1/NEMOGCM/TOOLS/SIREN/src/grid_zgr.f90 @ 7153

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

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

Download in other formats: