source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6986

Last change on this file since 6986 was 6986, checked in by acc, 5 years ago

Branch dev_r6393_NOC_WAD. Introduced some WAD test cases, tidied and corrected WAD code

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