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

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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domzgr.f90 @ 6951

Last change on this file since 6951 was 6951, checked in by flavoni, 8 years ago

merge simplif-2 branches (TOOLS and NEMO); update DOMAINcfg TOOL: create domain_cfg.nc files to be used in new version of NEMO, SIMPLIF-2 branch

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