New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.F90 in NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

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

Last change on this file since 11602 was 11602, checked in by mathiot, 5 years ago

ENHANCE-03_domcfg: remove gde3w_0 as outputed in domcfg.nc (ticket #2143)

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