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

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

ENHANCE-03_domcfg : add management of closed seas in domain cfg by flood filling and lat/lon seed instead of i/j box definition (ticket #2143)

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