source: branches/2017/dev_merge_2017/NEMOGCM/TOOLS/DOMAINcfg/src/domzgr.f90 @ 9033

Last change on this file since 9033 was 9033, checked in by timgraham, 3 years ago

Commit final files from merge of NEMOGCM and some fixes for waves (taooc renamed tauwoc)

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