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 @ 7025

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

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

Download in other formats: