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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domzgr.f90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/domzgr.f90 @ 10490

Last change on this file since 10490 was 10490, checked in by mathiot, 6 years ago

change the deepest e3w_0 inside the ice shelf => correct depth when using e3_to_depth

File size: 137.0 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 (./LICENSE)
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     (./LICENSE)
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               gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik)
1830               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
1831               e3w_0  (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2)
1832               gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1)
1833            ENDIF
1834         END DO
1835      END DO
1836   
1837      it = 0 
1838      DO jj = 1, jpj 
1839         DO ji = 1, jpi 
1840            ik = misfdep(ji,jj) 
1841            IF( ik > 1 ) THEN               ! ice shelf point only
1842               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
1843               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
1844            ! test
1845               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
1846               IF( zdiff <= 0. .AND. lwp ) THEN 
1847                  it = it + 1 
1848                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
1849                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
1850                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
1851                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
1852               ENDIF
1853            ENDIF
1854         END DO
1855      END DO
1856
1857      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
1858      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
1859      !
1860      IF( nn_timing == 1 )   CALL timing_stop('zgr_isf')
1861      !     
1862   END SUBROUTINE zgr_isf
1863
1864
1865   SUBROUTINE zgr_sco
1866      !!----------------------------------------------------------------------
1867      !!                  ***  ROUTINE zgr_sco  ***
1868      !!                     
1869      !! ** Purpose :   define the s-coordinate system
1870      !!
1871      !! ** Method  :   s-coordinate
1872      !!         The depth of model levels is defined as the product of an
1873      !!      analytical function by the local bathymetry, while the vertical
1874      !!      scale factors are defined as the product of the first derivative
1875      !!      of the analytical function by the bathymetry.
1876      !!      (this solution save memory as depth and scale factors are not
1877      !!      3d fields)
1878      !!          - Read bathymetry (in meters) at t-point and compute the
1879      !!         bathymetry at u-, v-, and f-points.
1880      !!            hbatu = mi( hbatt )
1881      !!            hbatv = mj( hbatt )
1882      !!            hbatf = mi( mj( hbatt ) )
1883      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1884      !!         function and its derivative given as function.
1885      !!            z_gsigt(k) = fssig (k    )
1886      !!            z_gsigw(k) = fssig (k-0.5)
1887      !!            z_esigt(k) = fsdsig(k    )
1888      !!            z_esigw(k) = fsdsig(k-0.5)
1889      !!      Three options for stretching are give, and they can be modified
1890      !!      following the users requirements. Nevertheless, the output as
1891      !!      well as the way to compute the model levels and scale factors
1892      !!      must be respected in order to insure second order accuracy
1893      !!      schemes.
1894      !!
1895      !!      The three methods for stretching available are:
1896      !!
1897      !!           s_sh94 (Song and Haidvogel 1994)
1898      !!                a sinh/tanh function that allows sigma and stretched sigma
1899      !!
1900      !!           s_sf12 (Siddorn and Furner 2012?)
1901      !!                allows the maintenance of fixed surface and or
1902      !!                bottom cell resolutions (cf. geopotential coordinates)
1903      !!                within an analytically derived stretched S-coordinate framework.
1904      !!
1905      !!          s_tanh  (Madec et al 1996)
1906      !!                a cosh/tanh function that gives stretched coordinates       
1907      !!
1908      !!----------------------------------------------------------------------
1909      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
1910      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
1911      INTEGER  ::   ios                      ! Local integer output status for namelist read
1912      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
1913      REAL(wp) ::   zrfact
1914      !
1915      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1916      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1917      !!
1918      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1919         &                rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1920     !!----------------------------------------------------------------------
1921      !
1922      IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
1923      !
1924      CALL wrk_alloc( jpi,jpj,   zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1925      !
1926      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1927      READ  ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1928901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1929
1930      REWIND( numnam_cfg )              ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1931      READ  ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1932902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1933      IF(lwm) WRITE ( numond, namzgr_sco )
1934
1935      IF(lwp) THEN                           ! control print
1936         WRITE(numout,*)
1937         WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1938         WRITE(numout,*) '~~~~~~~~~~~'
1939         WRITE(numout,*) '   Namelist namzgr_sco'
1940         WRITE(numout,*) '     stretching coeffs '
1941         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
1942         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
1943         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
1944         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
1945         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
1946         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
1947         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
1948         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
1949         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
1950         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
1951         WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
1952         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
1953         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
1954         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
1955         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
1956         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
1957         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
1958         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
1959      ENDIF
1960
1961      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
1962      hifu(:,:) = rn_sbot_min
1963      hifv(:,:) = rn_sbot_min
1964      hiff(:,:) = rn_sbot_min
1965
1966      !                                        ! set maximum ocean depth
1967      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1968
1969         DO jj = 1, jpj
1970            DO ji = 1, jpi
1971              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1972            END DO
1973         END DO
1974      !                                        ! =============================
1975      !                                        ! Define the envelop bathymetry   (hbatt)
1976      !                                        ! =============================
1977      ! use r-value to create hybrid coordinates
1978      zenv(:,:) = bathy(:,:)
1979      !
1980      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1981         DO jj = 1, jpj
1982            DO ji = 1, jpi
1983               IF( bathy(ji,jj) == 0._wp ) THEN
1984                  iip1 = MIN( ji+1, jpi )
1985                  ijp1 = MIN( jj+1, jpj )
1986                  iim1 = MAX( ji-1, 1 )
1987                  ijm1 = MAX( jj-1, 1 )
1988!!gm BUG fix see ticket #1617
1989                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1990                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1991                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1992                     &    zenv(ji,jj) = rn_sbot_min
1993!!gm
1994!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1995!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1996!!gm               zenv(ji,jj) = rn_sbot_min
1997!!gm             ENDIF
1998!!gm end
1999              ENDIF
2000            END DO
2001         END DO
2002
2003      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
2004      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2005      !
2006      ! smooth the bathymetry (if required)
2007      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
2008      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
2009      !
2010      jl = 0
2011      zrmax = 1._wp
2012      !   
2013      !     
2014      ! set scaling factor used in reducing vertical gradients
2015      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
2016      !
2017      ! initialise temporary evelope depth arrays
2018      ztmpi1(:,:) = zenv(:,:)
2019      ztmpi2(:,:) = zenv(:,:)
2020      ztmpj1(:,:) = zenv(:,:)
2021      ztmpj2(:,:) = zenv(:,:)
2022      !
2023      ! initialise temporary r-value arrays
2024      zri(:,:) = 1._wp
2025      zrj(:,:) = 1._wp
2026      !                                                            ! ================ !
2027      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
2028         !                                                         ! ================ !
2029         jl = jl + 1
2030         zrmax = 0._wp
2031         ! we set zrmax from previous r-values (zri and zrj) first
2032         ! if set after current r-value calculation (as previously)
2033         ! we could exit DO WHILE prematurely before checking r-value
2034         ! of current zenv
2035         DO jj = 1, nlcj
2036            DO ji = 1, nlci
2037               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2038            END DO
2039         END DO
2040         zri(:,:) = 0._wp
2041         zrj(:,:) = 0._wp
2042         DO jj = 1, nlcj
2043            DO ji = 1, nlci
2044               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
2045               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
2046               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2047                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
2048               END IF
2049               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2050                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
2051               END IF
2052               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
2053               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
2054               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
2055               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
2056            END DO
2057         END DO
2058         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
2059         !
2060         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
2061         !
2062         DO jj = 1, nlcj
2063            DO ji = 1, nlci
2064               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2065            END DO
2066         END DO
2067         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
2068         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2069         !                                                  ! ================ !
2070      END DO                                                !     End loop     !
2071      !                                                     ! ================ !
2072      DO jj = 1, jpj
2073         DO ji = 1, jpi
2074            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2075         END DO
2076      END DO
2077      !
2078      ! Envelope bathymetry saved in hbatt
2079      hbatt(:,:) = zenv(:,:) 
2080      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2081         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2082         DO jj = 1, jpj
2083            DO ji = 1, jpi
2084               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2085               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2086            END DO
2087         END DO
2088      ENDIF
2089      !
2090      !                                        ! ==============================
2091      !                                        !   hbatu, hbatv, hbatf fields
2092      !                                        ! ==============================
2093      IF(lwp) THEN
2094         WRITE(numout,*)
2095           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2096      ENDIF
2097      hbatu(:,:) = rn_sbot_min
2098      hbatv(:,:) = rn_sbot_min
2099      hbatf(:,:) = rn_sbot_min
2100      DO jj = 1, jpjm1
2101        DO ji = 1, jpim1   ! NO vector opt.
2102           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2103           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2104           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2105              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2106        END DO
2107      END DO
2108
2109      !
2110      ! Apply lateral boundary condition
2111!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2112      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2113      DO jj = 1, jpj
2114         DO ji = 1, jpi
2115            IF( hbatu(ji,jj) == 0._wp ) THEN
2116               !No worries about the following line when ln_wd == .true.
2117               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2118               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2119            ENDIF
2120         END DO
2121      END DO
2122      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2123      DO jj = 1, jpj
2124         DO ji = 1, jpi
2125            IF( hbatv(ji,jj) == 0._wp ) THEN
2126               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2127               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2128            ENDIF
2129         END DO
2130      END DO
2131      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2132      DO jj = 1, jpj
2133         DO ji = 1, jpi
2134            IF( hbatf(ji,jj) == 0._wp ) THEN
2135               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2136               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2137            ENDIF
2138         END DO
2139      END DO
2140
2141!!bug:  key_helsinki a verifer
2142        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2143        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2144        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2145        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2146
2147      IF( nprint == 1 .AND. lwp )   THEN
2148         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2149            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2150         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2151            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2152         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2153            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2154         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2155            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2156      ENDIF
2157!! helsinki
2158
2159      !                                            ! =======================
2160      !                                            !   s-ccordinate fields     (gdep., e3.)
2161      !                                            ! =======================
2162      !
2163      ! non-dimensional "sigma" for model level depth at w- and t-levels
2164
2165
2166!========================================================================
2167! Song and Haidvogel  1994 (ln_s_sh94=T)
2168! Siddorn and Furner 2012 (ln_sf12=T)
2169! or  tanh function       (both false)                   
2170!========================================================================
2171      IF      ( ln_s_sh94 ) THEN
2172                           CALL s_sh94()
2173      ELSE IF ( ln_s_sf12 ) THEN
2174                           CALL s_sf12()
2175      ELSE                 
2176                           CALL s_tanh()
2177      ENDIF
2178
2179      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2180      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2181      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2182      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2183      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2184      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2185      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2186      !
2187        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
2188        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
2189        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
2190        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
2191        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
2192        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
2193        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
2194
2195
2196!!gm   I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays)
2197!!gm   and only that !!!!!
2198!!gm   THIS should be removed from here !
2199      gdept_n(:,:,:) = gdept_0(:,:,:)
2200      gdepw_n(:,:,:) = gdepw_0(:,:,:)
2201      gde3w_n(:,:,:) = gde3w_0(:,:,:)
2202      e3t_n  (:,:,:) = e3t_0  (:,:,:)
2203      e3u_n  (:,:,:) = e3u_0  (:,:,:)
2204      e3v_n  (:,:,:) = e3v_0  (:,:,:)
2205      e3f_n  (:,:,:) = e3f_0  (:,:,:)
2206      e3w_n  (:,:,:) = e3w_0  (:,:,:)
2207      e3uw_n (:,:,:) = e3uw_0 (:,:,:)
2208      e3vw_n (:,:,:) = e3vw_0 (:,:,:)
2209!!gm and obviously in the following, use the _0 arrays until the end of this subroutine
2210!! gm end
2211!!
2212      ! HYBRID :
2213      DO jj = 1, jpj
2214         DO ji = 1, jpi
2215            DO jk = 1, jpkm1
2216               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2217            END DO
2218         END DO
2219      END DO
2220      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2221         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2222
2223      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2224         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2225         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2226            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gde3w_0(:,:,:) )
2227         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
2228            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
2229            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
2230            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2231
2232         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2233            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gde3w_0(:,:,:) )
2234         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
2235            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
2236            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
2237            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2238      ENDIF
2239      !  END DO
2240      IF(lwp) THEN                                  ! selected vertical profiles
2241         WRITE(numout,*)
2242         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2243         WRITE(numout,*) ' ~~~~~~  --------------------'
2244         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2245         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2246            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2247         DO jj = mj0(20), mj1(20)
2248            DO ji = mi0(20), mi1(20)
2249               WRITE(numout,*)
2250               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2251               WRITE(numout,*) ' ~~~~~~  --------------------'
2252               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2253               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2254                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2255            END DO
2256         END DO
2257         DO jj = mj0(74), mj1(74)
2258            DO ji = mi0(100), mi1(100)
2259               WRITE(numout,*)
2260               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2261               WRITE(numout,*) ' ~~~~~~  --------------------'
2262               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2263               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2264                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2265            END DO
2266         END DO
2267      ENDIF
2268      !
2269      !================================================================================
2270      ! check the coordinate makes sense
2271      !================================================================================
2272      DO ji = 1, jpi
2273         DO jj = 1, jpj
2274            !
2275            IF( hbatt(ji,jj) > 0._wp) THEN
2276               DO jk = 1, mbathy(ji,jj)
2277                 ! check coordinate is monotonically increasing
2278                 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN
2279                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2280                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2281                    WRITE(numout,*) 'e3w',e3w_n(ji,jj,:)
2282                    WRITE(numout,*) 'e3t',e3t_n(ji,jj,:)
2283                    CALL ctl_stop( ctmp1 )
2284                 ENDIF
2285                 ! and check it has never gone negative
2286                 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN
2287                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2288                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2289                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
2290                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2291                    CALL ctl_stop( ctmp1 )
2292                 ENDIF
2293                 ! and check it never exceeds the total depth
2294                 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2295                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2296                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2297                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
2298                    CALL ctl_stop( ctmp1 )
2299                 ENDIF
2300               END DO
2301               !
2302               DO jk = 1, mbathy(ji,jj)-1
2303                 ! and check it never exceeds the total depth
2304                IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2305                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2306                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2307                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2308                    CALL ctl_stop( ctmp1 )
2309                 ENDIF
2310               END DO
2311            ENDIF
2312         END DO
2313      END DO
2314      !
2315      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2316      !
2317      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2318      !
2319   END SUBROUTINE zgr_sco
2320
2321
2322   SUBROUTINE s_sh94()
2323      !!----------------------------------------------------------------------
2324      !!                  ***  ROUTINE s_sh94  ***
2325      !!                     
2326      !! ** Purpose :   stretch the s-coordinate system
2327      !!
2328      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2329      !!                mixed S/sigma coordinate
2330      !!
2331      !! Reference : Song and Haidvogel 1994.
2332      !!----------------------------------------------------------------------
2333      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2334      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2335      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
2336      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
2337      !
2338      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2339      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2340      !!----------------------------------------------------------------------
2341
2342      CALL wrk_alloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2343      CALL wrk_alloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2344
2345      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2346      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2347      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2348      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2349      !
2350      DO ji = 1, jpi
2351         DO jj = 1, jpj
2352            !
2353            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2354               DO jk = 1, jpk
2355                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2356                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2357               END DO
2358            ELSE ! shallow water, uniform sigma
2359               DO jk = 1, jpk
2360                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2361                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2362                  END DO
2363            ENDIF
2364            !
2365            DO jk = 1, jpkm1
2366               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2367               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2368            END DO
2369            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2370            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2371            !
2372            ! Coefficients for vertical depth as the sum of e3w scale factors
2373            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2374            DO jk = 2, jpk
2375               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2376            END DO
2377            !
2378            DO jk = 1, jpk
2379               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2380               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2381               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2382               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2383               gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2384            END DO
2385           !
2386         END DO   ! for all jj's
2387      END DO    ! for all ji's
2388
2389      DO ji = 1, jpim1
2390         DO jj = 1, jpjm1
2391            ! extended for Wetting/Drying case
2392            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
2393            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
2394            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
2395            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
2396            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
2397            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
2398                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
2399            DO jk = 1, jpk
2400                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
2401                        &              / ztmpu
2402                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
2403                        &              / ztmpu
2404
2405                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
2406                        &              / ztmpv
2407                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
2408                        &              / ztmpv
2409
2410                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
2411                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
2412                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
2413                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
2414
2415               !
2416               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2417               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2418               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2419               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2420               !
2421               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2422               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2423               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2424            END DO
2425        END DO
2426      END DO
2427      !
2428      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2429      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2430      !
2431   END SUBROUTINE s_sh94
2432
2433
2434   SUBROUTINE s_sf12
2435      !!----------------------------------------------------------------------
2436      !!                  ***  ROUTINE s_sf12 ***
2437      !!                     
2438      !! ** Purpose :   stretch the s-coordinate system
2439      !!
2440      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2441      !!                mixed S/sigma/Z coordinate
2442      !!
2443      !!                This method allows the maintenance of fixed surface and or
2444      !!                bottom cell resolutions (cf. geopotential coordinates)
2445      !!                within an analytically derived stretched S-coordinate framework.
2446      !!
2447      !!
2448      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2449      !!----------------------------------------------------------------------
2450      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2451      REAL(wp) ::   zsmth               ! smoothing around critical depth
2452      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2453      REAL(wp) ::   ztmpu, ztmpv, ztmpf
2454      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
2455      !
2456      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2457      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2458      !!----------------------------------------------------------------------
2459      !
2460      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2461      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2462
2463      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2464      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2465      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2466      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2467
2468      DO ji = 1, jpi
2469         DO jj = 1, jpj
2470
2471          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2472             
2473              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2474                                                     ! could be changed by users but care must be taken to do so carefully
2475              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2476           
2477              zzs = rn_zs / hbatt(ji,jj) 
2478             
2479              IF (rn_efold /= 0.0_wp) THEN
2480                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2481              ELSE
2482                zsmth = 1.0_wp 
2483              ENDIF
2484               
2485              DO jk = 1, jpk
2486                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2487                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2488              ENDDO
2489              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2490              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2491 
2492          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2493
2494            DO jk = 1, jpk
2495              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2496              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2497            END DO
2498
2499          ELSE  ! shallow water, z coordinates
2500
2501            DO jk = 1, jpk
2502              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2503              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2504            END DO
2505
2506          ENDIF
2507
2508          DO jk = 1, jpkm1
2509             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2510             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2511          END DO
2512          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2513          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2514
2515          ! Coefficients for vertical depth as the sum of e3w scale factors
2516          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2517          DO jk = 2, jpk
2518             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2519          END DO
2520
2521          DO jk = 1, jpk
2522             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2523             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2524             gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2525          END DO
2526
2527        ENDDO   ! for all jj's
2528      ENDDO    ! for all ji's
2529
2530      DO ji=1,jpi-1
2531        DO jj=1,jpj-1
2532
2533           ! extend to suit for Wetting/Drying case
2534           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
2535           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
2536           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
2537           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
2538           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
2539           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
2540                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
2541           DO jk = 1, jpk
2542                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
2543                       &              / ztmpu
2544                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
2545                       &              / ztmpu
2546                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
2547                       &              / ztmpv
2548                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
2549                       &              / ztmpv
2550
2551                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
2552                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
2553                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
2554                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
2555
2556             ! Code prior to wetting and drying option (for reference)
2557             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2558             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
2559             !
2560             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2561             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
2562             !
2563             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2564             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
2565             !
2566             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2567             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
2568             !
2569             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
2570             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
2571             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
2572             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
2573             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2574
2575             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2576             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2577             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2578             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2579             !
2580             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2581             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2582             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2583          END DO
2584
2585        ENDDO
2586      ENDDO
2587      !
2588      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2589      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2590      CALL lbc_lnk(e3w_0 ,'T',1.)
2591      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2592      !
2593      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2594      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2595      !
2596   END SUBROUTINE s_sf12
2597
2598
2599   SUBROUTINE s_tanh()
2600      !!----------------------------------------------------------------------
2601      !!                  ***  ROUTINE s_tanh***
2602      !!                     
2603      !! ** Purpose :   stretch the s-coordinate system
2604      !!
2605      !! ** Method  :   s-coordinate stretch
2606      !!
2607      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2608      !!----------------------------------------------------------------------
2609      INTEGER  ::   ji, jj, jk       ! dummy loop argument
2610      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2611      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2612      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2613      !!----------------------------------------------------------------------
2614
2615      CALL wrk_alloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
2616      CALL wrk_alloc( jpk,   z_esigt, z_esigw )
2617
2618      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2619      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2620
2621      DO jk = 1, jpk
2622        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2623        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2624      END DO
2625      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2626      !
2627      ! Coefficients for vertical scale factors at w-, t- levels
2628!!gm bug :  define it from analytical function, not like juste bellow....
2629!!gm        or betteroffer the 2 possibilities....
2630      DO jk = 1, jpkm1
2631         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2632         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2633      END DO
2634      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2635      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2636      !
2637      ! Coefficients for vertical depth as the sum of e3w scale factors
2638      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2639      DO jk = 2, jpk
2640         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2641      END DO
2642!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2643      DO jk = 1, jpk
2644         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2645         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2646         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2647         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2648         gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2649      END DO
2650!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2651      DO jj = 1, jpj
2652         DO ji = 1, jpi
2653            DO jk = 1, jpk
2654              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2655              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2656              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2657              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2658              !
2659              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2660              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2661              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2662            END DO
2663         END DO
2664      END DO
2665      !
2666      CALL wrk_dealloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
2667      CALL wrk_dealloc( jpk,   z_esigt, z_esigw          )
2668      !
2669   END SUBROUTINE s_tanh
2670
2671
2672   FUNCTION fssig( pk ) RESULT( pf )
2673      !!----------------------------------------------------------------------
2674      !!                 ***  ROUTINE fssig ***
2675      !!       
2676      !! ** Purpose :   provide the analytical function in s-coordinate
2677      !!         
2678      !! ** Method  :   the function provide the non-dimensional position of
2679      !!                T and W (i.e. between 0 and 1)
2680      !!                T-points at integer values (between 1 and jpk)
2681      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2682      !!----------------------------------------------------------------------
2683      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2684      REAL(wp)             ::   pf   ! sigma value
2685      !!----------------------------------------------------------------------
2686      !
2687      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2688         &     - TANH( rn_thetb * rn_theta                                )  )   &
2689         & * (   COSH( rn_theta                           )                      &
2690         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2691         & / ( 2._wp * SINH( rn_theta ) )
2692      !
2693   END FUNCTION fssig
2694
2695
2696   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2697      !!----------------------------------------------------------------------
2698      !!                 ***  ROUTINE fssig1 ***
2699      !!
2700      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2701      !!
2702      !! ** Method  :   the function provides the non-dimensional position of
2703      !!                T and W (i.e. between 0 and 1)
2704      !!                T-points at integer values (between 1 and jpk)
2705      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2706      !!----------------------------------------------------------------------
2707      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2708      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2709      REAL(wp)             ::   pf1   ! sigma value
2710      !!----------------------------------------------------------------------
2711      !
2712      IF ( rn_theta == 0 ) then      ! uniform sigma
2713         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2714      ELSE                        ! stretched sigma
2715         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2716            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2717            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2718      ENDIF
2719      !
2720   END FUNCTION fssig1
2721
2722
2723   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2724      !!----------------------------------------------------------------------
2725      !!                 ***  ROUTINE fgamma  ***
2726      !!
2727      !! ** Purpose :   provide analytical function for the s-coordinate
2728      !!
2729      !! ** Method  :   the function provides the non-dimensional position of
2730      !!                T and W (i.e. between 0 and 1)
2731      !!                T-points at integer values (between 1 and jpk)
2732      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2733      !!
2734      !!                This method allows the maintenance of fixed surface and or
2735      !!                bottom cell resolutions (cf. geopotential coordinates)
2736      !!                within an analytically derived stretched S-coordinate framework.
2737      !!
2738      !! Reference  :   Siddorn and Furner, in prep
2739      !!----------------------------------------------------------------------
2740      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2741      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2742      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
2743      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
2744      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
2745      !
2746      INTEGER  ::   jk             ! dummy loop index
2747      REAL(wp) ::   za1,za2,za3    ! local scalar
2748      REAL(wp) ::   zn1,zn2        !   -      -
2749      REAL(wp) ::   za,zb,zx       !   -      -
2750      !!----------------------------------------------------------------------
2751      !
2752      zn1  =  1._wp / REAL( jpkm1, wp )
2753      zn2  =  1._wp -  zn1
2754      !
2755      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2756      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2757      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2758      !
2759      za = pzb - za3*(pzs-za1)-za2
2760      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2761      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2762      zx = 1.0_wp-za/2.0_wp-zb
2763      !
2764      DO jk = 1, jpk
2765         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2766            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2767            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2768        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2769      END DO
2770      !
2771   END FUNCTION fgamma
2772
2773   !!======================================================================
2774END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.