source: utils/tools/SIREN/src/grid_zgr.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months ago

update nemo trunk

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