source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90 @ 11129

Last change on this file since 11129 was 11129, checked in by mathiot, 15 months ago

simplification of domcfg (rm all var_n and var_b as it is not needed) (ticket #2143)

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