source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6152

Last change on this file since 6152 was 6152, checked in by acc, 6 years ago

Add wetting and drying option from dev_r5803_NOC_WAD branch. Logically isolated code changes in domvvl.F90, domzgr.F90, dynhpg.F90, dynspg_ts.F90, sshwzv.F90 and nemogcm.F90. New module wet_dry.F90 in DYN. Fully SETTE tested with code deactivated (ln_wad=.false.). No test case yet available to justify activating option (still under development)

  • Property svn:keywords set to Id
File size: 139.6 KB
Line 
1MODULE domzgr
2   !!==============================================================================
3   !!                       ***  MODULE domzgr   ***
4   !! Ocean domain : definition of the vertical coordinate system
5   !!==============================================================================
6   !! History :  OPA  ! 1995-12  (G. Madec)  Original code : s vertical coordinate
7   !!                 ! 1997-07  (G. Madec)  lbc_lnk call
8   !!                 ! 1997-04  (J.-O. Beismann)
9   !!            8.5  ! 2002-09  (A. Bozec, G. Madec)  F90: Free form and module
10   !!             -   ! 2002-09  (A. de Miranda)  rigid-lid + islands
11   !!  NEMO      1.0  ! 2003-08  (G. Madec)  F90: Free form and module
12   !!             -   ! 2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
13   !!            2.0  ! 2006-04  (R. Benshila, G. Madec)  add zgr_zco
14   !!            3.0  ! 2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
15   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
16   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
18   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case 
19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
20   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dom_zgr          : defined the ocean vertical coordinate system
25   !!       zgr_bat      : bathymetry fields (levels and meters)
26   !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
27   !!       zgr_bat_ctl  : check the bathymetry files
28   !!       zgr_bot_level: deepest ocean level for t-, u, and v-points
29   !!       zgr_z        : reference z-coordinate
30   !!       zgr_zco      : z-coordinate
31   !!       zgr_zps      : z-coordinate with partial steps
32   !!       zgr_sco      : s-coordinate
33   !!       fssig        : tanh stretch function
34   !!       fssig1       : Song and Haidvogel 1994 stretch function
35   !!       fgamma       : Siddorn and Furner 2012 stretching function
36   !!---------------------------------------------------------------------
37   USE oce               ! ocean variables
38   USE dom_oce           ! ocean domain
39   USE wet_dry           ! wetting and drying
40   USE closea            ! closed seas
41   USE c1d               ! 1D vertical configuration
42   !
43   USE in_out_manager    ! I/O manager
44   USE iom               ! I/O library
45   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
46   USE lib_mpp           ! distributed memory computing library
47   USE wrk_nemo          ! Memory allocation
48   USE timing            ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   dom_zgr        ! called by dom_init.F90
54
55   !                              !!* Namelist namzgr_sco *
56   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
57   LOGICAL  ::   ln_s_sf12         ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
58   !
59   REAL(wp) ::   rn_sbot_min       ! minimum depth of s-bottom surface (>0) (m)
60   REAL(wp) ::   rn_sbot_max       ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
61   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1)
62   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates
63   ! Song and Haidvogel 1994 stretching parameters
64   REAL(wp) ::   rn_theta          ! surface control parameter (0<=rn_theta<=20)
65   REAL(wp) ::   rn_thetb          ! bottom control parameter  (0<=rn_thetb<= 1)
66   REAL(wp) ::   rn_bb             ! stretching parameter
67   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
68   ! Siddorn and Furner stretching parameters
69   LOGICAL  ::   ln_sigcrit        ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
70   REAL(wp) ::   rn_alpha          ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
71   REAL(wp) ::   rn_efold          !  efold length scale for transition to stretched coord
72   REAL(wp) ::   rn_zs             !  depth of surface grid box
73                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
74   REAL(wp) ::   rn_zb_a           !  bathymetry scaling factor for calculating Zb
75   REAL(wp) ::   rn_zb_b           !  offset for calculating Zb
76
77  !! * Substitutions
78#  include "vectopt_loop_substitute.h90"
79   !!----------------------------------------------------------------------
80   !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
81   !! $Id$
82   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
83   !!----------------------------------------------------------------------
84CONTAINS       
85
86   SUBROUTINE dom_zgr
87      !!----------------------------------------------------------------------
88      !!                ***  ROUTINE dom_zgr  ***
89      !!                   
90      !! ** Purpose :   set the depth of model levels and the resulting
91      !!              vertical scale factors.
92      !!
93      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
94      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
95      !!              - vertical coordinate (gdep., e3.) depending on the
96      !!                coordinate chosen :
97      !!                   ln_zco=T   z-coordinate   
98      !!                   ln_zps=T   z-coordinate with partial steps
99      !!                   ln_zco=T   s-coordinate
100      !!
101      !! ** Action  :   define gdep., e3., mbathy and bathy
102      !!----------------------------------------------------------------------
103      INTEGER ::   ioptio, ibat   ! local integer
104      INTEGER ::   ios
105      !
106      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
107      !!----------------------------------------------------------------------
108      !
109      IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
110      !
111      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
112      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
113901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
114
115      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
116      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
117902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
118      IF(lwm) WRITE ( numond, namzgr )
119
120      IF(lwp) THEN                     ! Control print
121         WRITE(numout,*)
122         WRITE(numout,*) 'dom_zgr : vertical coordinate'
123         WRITE(numout,*) '~~~~~~~'
124         WRITE(numout,*) '   Namelist namzgr : set vertical coordinate'
125         WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
126         WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
127         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco
128         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav
129         WRITE(numout,*) '      linear free surface            ln_linssh = ', ln_linssh
130      ENDIF
131
132      IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time'
133
134      ioptio = 0                       ! Check Vertical coordinate options
135      IF( ln_zco      )   ioptio = ioptio + 1
136      IF( ln_zps      )   ioptio = ioptio + 1
137      IF( ln_sco      )   ioptio = ioptio + 1
138      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
139      !
140      ! Build the vertical coordinate system
141      ! ------------------------------------
142                          CALL zgr_z            ! Reference z-coordinate system (always called)
143                          CALL zgr_bat          ! Bathymetry fields (levels and meters)
144      IF( lk_c1d      )   CALL lbc_lnk( bathy , 'T', 1._wp )   ! 1D config.: same bathy value over the 3x3 domain
145      IF( ln_zco      )   CALL zgr_zco          ! z-coordinate
146      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate
147      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate
148      !
149      ! final adjustment of mbathy & check
150      ! -----------------------------------
151      IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain
152      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points
153                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points
154                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points
155      !
156      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain
157         ibat = mbathy(2,2)
158         mbathy(:,:) = ibat
159      END IF
160      !
161      IF( nprint == 1 .AND. lwp )   THEN
162         WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
163         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
164            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
165         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  &
166            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  &
167            &                          ' uw', MINVAL(  e3uw_0(:,:,:) ), ' vw', MINVAL(  e3vw_0(:,:,:)),   &
168            &                          ' w ', MINVAL(   e3w_0(:,:,:) )
169
170         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
171            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
172         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  &
173            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  &
174            &                          ' uw', MAXVAL(  e3uw_0(:,:,:) ), ' vw', MAXVAL(  e3vw_0(:,:,:) ),  &
175            &                          ' w ', MAXVAL(   e3w_0(:,:,:) )
176      ENDIF
177      !
178      IF( nn_timing == 1 )  CALL timing_stop('dom_zgr')
179      !
180   END SUBROUTINE dom_zgr
181
182
183   SUBROUTINE zgr_z
184      !!----------------------------------------------------------------------
185      !!                   ***  ROUTINE zgr_z  ***
186      !!                   
187      !! ** Purpose :   set the depth of model levels and the resulting
188      !!      vertical scale factors.
189      !!
190      !! ** Method  :   z-coordinate system (use in all type of coordinate)
191      !!        The depth of model levels is defined from an analytical
192      !!      function the derivative of which gives the scale factors.
193      !!        both depth and scale factors only depend on k (1d arrays).
194      !!              w-level: gdepw_1d  = gdep(k)
195      !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
196      !!              t-level: gdept_1d  = gdep(k+0.5)
197      !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
198      !!
199      !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
200      !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
201      !!
202      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
203      !!----------------------------------------------------------------------
204      INTEGER  ::   jk                     ! dummy loop indices
205      REAL(wp) ::   zt, zw                 ! temporary scalars
206      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
207      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
208      REAL(wp) ::   zrefdep                ! depth of the reference level (~10m)
209      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
210      !!----------------------------------------------------------------------
211      !
212      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
213      !
214      ! Set variables from parameters
215      ! ------------------------------
216       zkth = ppkth       ;   zacr = ppacr
217       zdzmin = ppdzmin   ;   zhmax = pphmax
218       zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
219
220      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
221      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
222      IF(   ppa1  == pp_to_be_computed  .AND.  &
223         &  ppa0  == pp_to_be_computed  .AND.  &
224         &  ppsur == pp_to_be_computed           ) THEN
225         !
226#if defined key_agrif
227         za1  = (  ppdzmin - pphmax / FLOAT(jpkdta-1)  )                                                   &
228            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * (  LOG( COSH( (jpkdta - ppkth) / ppacr) )&
229            &                                                      - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
230#else
231         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
232            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
233            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
234#endif
235         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
236         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
237      ELSE
238         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
239         za2 = ppa2                            ! optional (ldbletanh=T) double tanh parameter
240      ENDIF
241
242      IF(lwp) THEN                         ! Parameter print
243         WRITE(numout,*)
244         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
245         WRITE(numout,*) '    ~~~~~~~'
246         IF(  ppkth == 0._wp ) THEN             
247              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
248              WRITE(numout,*) '            Total depth    :', zhmax
249#if defined key_agrif
250              WRITE(numout,*) '            Layer thickness:', zhmax/(jpkdta-1)
251#else
252              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
253#endif
254         ELSE
255            IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
256               WRITE(numout,*) '         zsur, za0, za1 computed from '
257               WRITE(numout,*) '                 zdzmin = ', zdzmin
258               WRITE(numout,*) '                 zhmax  = ', zhmax
259            ENDIF
260            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
261            WRITE(numout,*) '                 zsur = ', zsur
262            WRITE(numout,*) '                 za0  = ', za0
263            WRITE(numout,*) '                 za1  = ', za1
264            WRITE(numout,*) '                 zkth = ', zkth
265            WRITE(numout,*) '                 zacr = ', zacr
266            IF( ldbletanh ) THEN
267               WRITE(numout,*) ' (Double tanh    za2  = ', za2
268               WRITE(numout,*) '  parameters)    zkth2= ', zkth2
269               WRITE(numout,*) '                 zacr2= ', zacr2
270            ENDIF
271         ENDIF
272      ENDIF
273
274
275      ! Reference z-coordinate (depth - scale factor at T- and W-points)
276      ! ======================
277      IF( ppkth == 0._wp ) THEN            !  uniform vertical grid
278#if defined key_agrif
279         za1 = zhmax / FLOAT(jpkdta-1) 
280#else
281         za1 = zhmax / FLOAT(jpk-1) 
282#endif
283         DO jk = 1, jpk
284            zw = FLOAT( jk )
285            zt = FLOAT( jk ) + 0.5_wp
286            gdepw_1d(jk) = ( zw - 1 ) * za1
287            gdept_1d(jk) = ( zt - 1 ) * za1
288            e3w_1d  (jk) =  za1
289            e3t_1d  (jk) =  za1
290         END DO
291      ELSE                                ! Madec & Imbard 1996 function
292         IF( .NOT. ldbletanh ) THEN
293            DO jk = 1, jpk
294               zw = REAL( jk , wp )
295               zt = REAL( jk , wp ) + 0.5_wp
296               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
297               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
298               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
299               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
300            END DO
301         ELSE
302            DO jk = 1, jpk
303               zw = FLOAT( jk )
304               zt = FLOAT( jk ) + 0.5_wp
305               ! Double tanh function
306               gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
307                  &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
308               gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
309                  &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
310               e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
311                  &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
312               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
313                  &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
314            END DO
315         ENDIF
316         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
317      ENDIF
318
319      IF ( ln_isfcav ) THEN
320! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
321! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
322         DO jk = 1, jpkm1
323            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
324         END DO
325         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
326
327         DO jk = 2, jpk
328            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
329         END DO
330         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
331      END IF
332
333!!gm BUG in s-coordinate this does not work!
334      ! deepest/shallowest W level Above/Below ~10m
335      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness)
336      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
337      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
338!!gm end bug
339
340      IF(lwp) THEN                        ! control print
341         WRITE(numout,*)
342         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
343         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
344         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
345      ENDIF
346      DO jk = 1, jpk                      ! control positivity
347         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
348         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
349      END DO
350      !
351      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
352      !
353   END SUBROUTINE zgr_z
354
355
356   SUBROUTINE zgr_bat
357      !!----------------------------------------------------------------------
358      !!                    ***  ROUTINE zgr_bat  ***
359      !!
360      !! ** Purpose :   set bathymetry both in levels and meters
361      !!
362      !! ** Method  :   read or define mbathy and bathy arrays
363      !!       * level bathymetry:
364      !!      The ocean basin geometry is given by a two-dimensional array,
365      !!      mbathy, which is defined as follow :
366      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
367      !!                              at t-point (ji,jj).
368      !!                            = 0  over the continental t-point.
369      !!      The array mbathy is checked to verified its consistency with
370      !!      model option. in particular:
371      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
372      !!                  along closed boundary.
373      !!            mbathy must be cyclic IF jperio=1.
374      !!            mbathy must be lower or equal to jpk-1.
375      !!            isolated ocean grid points are suppressed from mbathy
376      !!                  since they are only connected to remaining
377      !!                  ocean through vertical diffusion.
378      !!      ntopo=-1 :   rectangular channel or bassin with a bump
379      !!      ntopo= 0 :   flat rectangular channel or basin
380      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
381      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
382      !!
383      !! ** Action  : - mbathy: level bathymetry (in level index)
384      !!              - bathy : meter bathymetry (in meters)
385      !!----------------------------------------------------------------------
386      INTEGER  ::   ji, jj, jk            ! dummy loop indices
387      INTEGER  ::   inum                      ! temporary logical unit
388      INTEGER  ::   ierror                    ! error flag
389      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
390      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
391      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics
392      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
393      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data
394      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data
395      !!----------------------------------------------------------------------
396      !
397      IF( nn_timing == 1 )  CALL timing_start('zgr_bat')
398      !
399      IF(lwp) WRITE(numout,*)
400      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
401      IF(lwp) WRITE(numout,*) '    ~~~~~~~'
402      !                                               ! ================== !
403      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
404         !                                            ! ================== !
405         !                                            ! global domain level and meter bathymetry (idta,zdta)
406         !
407         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror )
408         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
409         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror )
410         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
411         !
412         IF( ntopo == 0 ) THEN                        ! flat basin
413            IF(lwp) WRITE(numout,*)
414            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
415            IF( rn_bathy > 0.01 ) THEN
416               IF(lwp) WRITE(numout,*) '         Depth = rn_bathy read in namelist'
417               zdta(:,:) = rn_bathy
418               IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
419                  idta(:,:) = jpkm1
420               ELSE                                                ! z-coordinate (zco or zps): step-like topography
421                  idta(:,:) = jpkm1
422                  DO jk = 1, jpkm1
423                     WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
424                  END DO
425               ENDIF
426            ELSE
427               IF(lwp) WRITE(numout,*) '         Depth = depthw(jpkm1)'
428               idta(:,:) = jpkm1                            ! before last level
429               zdta(:,:) = gdepw_1d(jpk)                     ! last w-point depth
430               h_oce     = gdepw_1d(jpk)
431            ENDIF
432         ELSE                                         ! bump centered in the basin
433            IF(lwp) WRITE(numout,*)
434            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
435            ii_bump = jpidta / 2                           ! i-index of the bump center
436            ij_bump = jpjdta / 2                           ! j-index of the bump center
437            r_bump  = 50000._wp                            ! bump radius (meters)       
438            h_bump  =  2700._wp                            ! bump height (meters)
439            h_oce   = gdepw_1d(jpk)                        ! background ocean depth (meters)
440            IF(lwp) WRITE(numout,*) '            bump characteristics: '
441            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
442            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
443            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
444            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
445            !                                       
446            DO jj = 1, jpjdta                              ! zdta :
447               DO ji = 1, jpidta
448                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
449                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
450                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
451               END DO
452            END DO
453            !                                              ! idta :
454            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
455               idta(:,:) = jpkm1
456            ELSE                                                ! z-coordinate (zco or zps): step-like topography
457               idta(:,:) = jpkm1
458               DO jk = 1, jpkm1
459                  WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) )   idta(:,:) = jk
460               END DO
461            ENDIF
462         ENDIF
463         !                                            ! set GLOBAL boundary conditions
464         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
465         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
466            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp
467            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp
468         ELSEIF( jperio == 2 ) THEN
469            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
470            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp
471            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp
472            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp
473         ELSE
474            ih = 0                                  ;      zh = 0._wp
475            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
476            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
477            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
478            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
479            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
480         ENDIF
481
482         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
483         mbathy(:,:) = 0                                   ! set to zero extra halo points
484         bathy (:,:) = 0._wp                               ! (require for mpp case)
485         DO jj = 1, nlcj                                   ! interior values
486            DO ji = 1, nlci
487               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
488               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
489            END DO
490         END DO
491         risfdep(:,:)=0.e0
492         misfdep(:,:)=1
493         !
494         DEALLOCATE( idta, zdta )
495         !
496         !                                            ! ================ !
497      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
498         !                                            ! ================ !
499         !
500         IF( ln_zco )   THEN                          ! zco : read level bathymetry
501            CALL iom_open ( 'bathy_level.nc', inum ) 
502            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
503            CALL iom_close( inum )
504            mbathy(:,:) = INT( bathy(:,:) )
505            !                                                ! =====================
506            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
507               !                                             ! =====================
508               !
509               ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open
510               ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
511               DO ji = mi0(ii0), mi1(ii1)
512                  DO jj = mj0(ij0), mj1(ij1)
513                     mbathy(ji,jj) = 15
514                  END DO
515               END DO
516               IF(lwp) WRITE(numout,*)
517               IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
518               !
519               ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
520               ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
521               DO ji = mi0(ii0), mi1(ii1)
522                  DO jj = mj0(ij0), mj1(ij1)
523                     mbathy(ji,jj) = 12
524                  END DO
525               END DO
526               IF(lwp) WRITE(numout,*)
527               IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
528               !
529            ENDIF
530            !
531         ENDIF
532         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
533            CALL iom_open ( 'bathy_meter.nc', inum ) 
534            IF ( ln_isfcav ) THEN
535               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
536            ELSE
537               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  )
538            END IF
539            CALL iom_close( inum )
540            !                                               
541            risfdep(:,:)=0._wp         
542            misfdep(:,:)=1             
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, fs_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, fs_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, fs_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      IF( .NOT.ln_wd ) THEN                     
1945         DO jj = 1, jpj
1946            DO ji = 1, jpi
1947              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1948            END DO
1949         END DO
1950      END IF
1951      !                                        ! =============================
1952      !                                        ! Define the envelop bathymetry   (hbatt)
1953      !                                        ! =============================
1954      ! use r-value to create hybrid coordinates
1955      zenv(:,:) = bathy(:,:)
1956      !
1957      IF( .NOT.ln_wd ) THEN   
1958      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1959         DO jj = 1, jpj
1960            DO ji = 1, jpi
1961               IF( bathy(ji,jj) == 0._wp ) THEN
1962                  iip1 = MIN( ji+1, jpi )
1963                  ijp1 = MIN( jj+1, jpj )
1964                  iim1 = MAX( ji-1, 1 )
1965                  ijm1 = MAX( jj-1, 1 )
1966!!gm BUG fix see ticket #1617
1967                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              &
1968                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              &
1969                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) &
1970                     &    zenv(ji,jj) = rn_sbot_min
1971!!gm
1972!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         &
1973!!gm                  &  bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1974!!gm               zenv(ji,jj) = rn_sbot_min
1975!!gm             ENDIF
1976!!gm end
1977              ENDIF
1978            END DO
1979         END DO
1980      END IF
1981
1982      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
1983      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
1984      !
1985      ! smooth the bathymetry (if required)
1986      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
1987      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
1988      !
1989      jl = 0
1990      zrmax = 1._wp
1991      !   
1992      !     
1993      ! set scaling factor used in reducing vertical gradients
1994      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1995      !
1996      ! initialise temporary evelope depth arrays
1997      ztmpi1(:,:) = zenv(:,:)
1998      ztmpi2(:,:) = zenv(:,:)
1999      ztmpj1(:,:) = zenv(:,:)
2000      ztmpj2(:,:) = zenv(:,:)
2001      !
2002      ! initialise temporary r-value arrays
2003      zri(:,:) = 1._wp
2004      zrj(:,:) = 1._wp
2005      !                                                            ! ================ !
2006      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
2007         !                                                         ! ================ !
2008         jl = jl + 1
2009         zrmax = 0._wp
2010         ! we set zrmax from previous r-values (zri and zrj) first
2011         ! if set after current r-value calculation (as previously)
2012         ! we could exit DO WHILE prematurely before checking r-value
2013         ! of current zenv
2014         DO jj = 1, nlcj
2015            DO ji = 1, nlci
2016               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
2017            END DO
2018         END DO
2019         zri(:,:) = 0._wp
2020         zrj(:,:) = 0._wp
2021         DO jj = 1, nlcj
2022            DO ji = 1, nlci
2023               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
2024               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
2025               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
2026                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
2027               END IF
2028               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
2029                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
2030               END IF
2031               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
2032               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
2033               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
2034               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
2035            END DO
2036         END DO
2037         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
2038         !
2039         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
2040         !
2041         DO jj = 1, nlcj
2042            DO ji = 1, nlci
2043               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
2044            END DO
2045         END DO
2046         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
2047         CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
2048         !                                                  ! ================ !
2049      END DO                                                !     End loop     !
2050      !                                                     ! ================ !
2051      DO jj = 1, jpj
2052         DO ji = 1, jpi
2053            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
2054         END DO
2055      END DO
2056      !
2057      ! Envelope bathymetry saved in hbatt
2058      hbatt(:,:) = zenv(:,:) 
2059      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2060         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2061         DO jj = 1, jpj
2062            DO ji = 1, jpi
2063               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
2064               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2065            END DO
2066         END DO
2067      ENDIF
2068      !
2069      !                                        ! ==============================
2070      !                                        !   hbatu, hbatv, hbatf fields
2071      !                                        ! ==============================
2072      IF(lwp) THEN
2073         WRITE(numout,*)
2074         IF( .NOT.ln_wd ) THEN
2075           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
2076         ELSE
2077           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min
2078           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld
2079         ENDIF
2080      ENDIF
2081      hbatu(:,:) = rn_sbot_min
2082      hbatv(:,:) = rn_sbot_min
2083      hbatf(:,:) = rn_sbot_min
2084      DO jj = 1, jpjm1
2085        DO ji = 1, jpim1   ! NO vector opt.
2086           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
2087           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
2088           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
2089              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
2090        END DO
2091      END DO
2092
2093      IF( ln_wd ) THEN               !avoid the zero depth on T- (U-,V-,F-) points
2094        DO jj = 1, jpj
2095          DO ji = 1, jpi
2096            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &
2097              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1
2098            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &
2099              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1
2100            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &
2101              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1
2102            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &
2103              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1
2104          END DO
2105        END DO
2106      END IF
2107      !
2108      ! Apply lateral boundary condition
2109!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
2110      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
2111      DO jj = 1, jpj
2112         DO ji = 1, jpi
2113            IF( hbatu(ji,jj) == 0._wp ) THEN
2114               !No worries about the following line when ln_wd == .true.
2115               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
2116               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
2117            ENDIF
2118         END DO
2119      END DO
2120      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
2121      DO jj = 1, jpj
2122         DO ji = 1, jpi
2123            IF( hbatv(ji,jj) == 0._wp ) THEN
2124               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
2125               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
2126            ENDIF
2127         END DO
2128      END DO
2129      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
2130      DO jj = 1, jpj
2131         DO ji = 1, jpi
2132            IF( hbatf(ji,jj) == 0._wp ) THEN
2133               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
2134               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
2135            ENDIF
2136         END DO
2137      END DO
2138
2139!!bug:  key_helsinki a verifer
2140      IF( .NOT.ln_wd ) THEN
2141        hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
2142        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
2143        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
2144        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
2145      END IF
2146
2147      IF( nprint == 1 .AND. lwp )   THEN
2148         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
2149            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
2150         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
2151            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
2152         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
2153            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
2154         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
2155            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
2156      ENDIF
2157!! helsinki
2158
2159      !                                            ! =======================
2160      !                                            !   s-ccordinate fields     (gdep., e3.)
2161      !                                            ! =======================
2162      !
2163      ! non-dimensional "sigma" for model level depth at w- and t-levels
2164
2165
2166!========================================================================
2167! Song and Haidvogel  1994 (ln_s_sh94=T)
2168! Siddorn and Furner 2012 (ln_sf12=T)
2169! or  tanh function       (both false)                   
2170!========================================================================
2171      IF      ( ln_s_sh94 ) THEN
2172                           CALL s_sh94()
2173      ELSE IF ( ln_s_sf12 ) THEN
2174                           CALL s_sf12()
2175      ELSE                 
2176                           CALL s_tanh()
2177      ENDIF
2178
2179      CALL lbc_lnk( e3t_0 , 'T', 1._wp )
2180      CALL lbc_lnk( e3u_0 , 'U', 1._wp )
2181      CALL lbc_lnk( e3v_0 , 'V', 1._wp )
2182      CALL lbc_lnk( e3f_0 , 'F', 1._wp )
2183      CALL lbc_lnk( e3w_0 , 'W', 1._wp )
2184      CALL lbc_lnk( e3uw_0, 'U', 1._wp )
2185      CALL lbc_lnk( e3vw_0, 'V', 1._wp )
2186      !
2187      IF( .NOT.ln_wd ) THEN
2188        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp
2189        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp
2190        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp
2191        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp
2192        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp
2193        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp
2194        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp
2195      END IF
2196
2197#if defined key_agrif
2198      IF( .NOT. Agrif_Root() ) THEN    ! Ensure meaningful vertical scale factors in ghost lines/columns
2199         IF( nbondi == -1 .OR. nbondi == 2 )   e3u_0(  1   ,  :   ,:) = e3u_0(  2   ,  :   ,:)
2200         IF( nbondi ==  1 .OR. nbondi == 2 )   e3u_0(nlci-1,  :   ,:) = e3u_0(nlci-2,  :   ,:)
2201         IF( nbondj == -1 .OR. nbondj == 2 )   e3v_0(  :   ,  1   ,:) = e3v_0(  :   ,  2   ,:)
2202         IF( nbondj ==  1 .OR. nbondj == 2 )   e3v_0(  :   ,nlcj-1,:) = e3v_0(  :   ,nlcj-2,:)
2203       ENDIF
2204#endif
2205
2206!!gm   I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays)
2207!!gm   and only that !!!!!
2208!!gm   THIS should be removed from here !
2209      gdept_n(:,:,:) = gdept_0(:,:,:)
2210      gdepw_n(:,:,:) = gdepw_0(:,:,:)
2211      gde3w_n(:,:,:) = gde3w_0(:,:,:)
2212      e3t_n  (:,:,:) = e3t_0  (:,:,:)
2213      e3u_n  (:,:,:) = e3u_0  (:,:,:)
2214      e3v_n  (:,:,:) = e3v_0  (:,:,:)
2215      e3f_n  (:,:,:) = e3f_0  (:,:,:)
2216      e3w_n  (:,:,:) = e3w_0  (:,:,:)
2217      e3uw_n (:,:,:) = e3uw_0 (:,:,:)
2218      e3vw_n (:,:,:) = e3vw_0 (:,:,:)
2219!!gm and obviously in the following, use the _0 arrays until the end of this subroutine
2220!! gm end
2221!!
2222      ! HYBRID :
2223      DO jj = 1, jpj
2224         DO ji = 1, jpi
2225            DO jk = 1, jpkm1
2226               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
2227            END DO
2228            IF( ln_wd ) THEN
2229              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
2230                mbathy(ji,jj) = 0
2231              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN
2232                mbathy(ji,jj) = 2
2233              ENDIF
2234            ELSE
2235              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0
2236            ENDIF
2237         END DO
2238      END DO
2239      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
2240         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
2241
2242      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
2243         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
2244         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
2245            &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gde3w_0(:,:,:) )
2246         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   &
2247            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   &
2248            &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0 (:,:,:) ),   &
2249            &                          ' w ', MINVAL( e3w_0  (:,:,:) )
2250
2251         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
2252            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gde3w_0(:,:,:) )
2253         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   &
2254            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   &
2255            &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0 (:,:,:) ),   &
2256            &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
2257      ENDIF
2258      !  END DO
2259      IF(lwp) THEN                                  ! selected vertical profiles
2260         WRITE(numout,*)
2261         WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
2262         WRITE(numout,*) ' ~~~~~~  --------------------'
2263         WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2264         WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
2265            &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
2266         DO jj = mj0(20), mj1(20)
2267            DO ji = mi0(20), mi1(20)
2268               WRITE(numout,*)
2269               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2270               WRITE(numout,*) ' ~~~~~~  --------------------'
2271               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2272               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2273                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2274            END DO
2275         END DO
2276         DO jj = mj0(74), mj1(74)
2277            DO ji = mi0(100), mi1(100)
2278               WRITE(numout,*)
2279               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
2280               WRITE(numout,*) ' ~~~~~~  --------------------'
2281               WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
2282               WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
2283                  &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
2284            END DO
2285         END DO
2286      ENDIF
2287      !
2288      !================================================================================
2289      ! check the coordinate makes sense
2290      !================================================================================
2291      DO ji = 1, jpi
2292         DO jj = 1, jpj
2293            !
2294            IF( hbatt(ji,jj) > 0._wp) THEN
2295               DO jk = 1, mbathy(ji,jj)
2296                 ! check coordinate is monotonically increasing
2297                 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN
2298                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2299                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
2300                    WRITE(numout,*) 'e3w',e3w_n(ji,jj,:)
2301                    WRITE(numout,*) 'e3t',e3t_n(ji,jj,:)
2302                    CALL ctl_stop( ctmp1 )
2303                 ENDIF
2304                 ! and check it has never gone negative
2305                 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN
2306                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
2307                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
2308                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
2309                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2310                    CALL ctl_stop( ctmp1 )
2311                 ENDIF
2312                 ! and check it never exceeds the total depth
2313                 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2314                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2315                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
2316                    WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:)
2317                    CALL ctl_stop( ctmp1 )
2318                 ENDIF
2319               END DO
2320               !
2321               DO jk = 1, mbathy(ji,jj)-1
2322                 ! and check it never exceeds the total depth
2323                IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN
2324                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2325                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
2326                    WRITE(numout,*) 'gdept',gdept_n(ji,jj,:)
2327                    CALL ctl_stop( ctmp1 )
2328                 ENDIF
2329               END DO
2330            ENDIF
2331         END DO
2332      END DO
2333      !
2334      CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
2335      !
2336      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
2337      !
2338   END SUBROUTINE zgr_sco
2339
2340
2341   SUBROUTINE s_sh94()
2342      !!----------------------------------------------------------------------
2343      !!                  ***  ROUTINE s_sh94  ***
2344      !!                     
2345      !! ** Purpose :   stretch the s-coordinate system
2346      !!
2347      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
2348      !!                mixed S/sigma coordinate
2349      !!
2350      !! Reference : Song and Haidvogel 1994.
2351      !!----------------------------------------------------------------------
2352      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2353      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2354      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf
2355      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
2356      !
2357      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2358      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2359      !!----------------------------------------------------------------------
2360
2361      CALL wrk_alloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2362      CALL wrk_alloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2363
2364      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2365      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2366      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2367      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2368      !
2369      DO ji = 1, jpi
2370         DO jj = 1, jpj
2371            !
2372            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
2373               DO jk = 1, jpk
2374                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
2375                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
2376               END DO
2377            ELSE ! shallow water, uniform sigma
2378               DO jk = 1, jpk
2379                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
2380                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2381                  END DO
2382            ENDIF
2383            !
2384            DO jk = 1, jpkm1
2385               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2386               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2387            END DO
2388            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
2389            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
2390            !
2391            ! Coefficients for vertical depth as the sum of e3w scale factors
2392            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
2393            DO jk = 2, jpk
2394               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2395            END DO
2396            !
2397            DO jk = 1, jpk
2398               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2399               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2400               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
2401               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
2402               gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
2403            END DO
2404           !
2405         END DO   ! for all jj's
2406      END DO    ! for all ji's
2407
2408      DO ji = 1, jpim1
2409         DO jj = 1, jpjm1
2410            ! extended for Wetting/Drying case
2411            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
2412            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
2413            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
2414            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
2415            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
2416            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
2417                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
2418            DO jk = 1, jpk
2419               IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
2420                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
2421                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
2422               ELSE
2423                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
2424                        &              / ztmpu
2425                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
2426                        &              / ztmpu
2427               END IF
2428
2429               IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
2430                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
2431                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
2432               ELSE
2433                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
2434                        &              / ztmpv
2435                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
2436                        &              / ztmpv
2437               END IF
2438
2439               IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
2440                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj  ,jk) + z_esigt3(ji+1,jj  ,jk)               &
2441                        &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
2442               ELSE
2443                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
2444                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
2445                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
2446                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
2447               END IF
2448
2449               !
2450               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2451               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2452               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2453               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2454               !
2455               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2456               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2457               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
2458            END DO
2459        END DO
2460      END DO
2461      !
2462      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2463      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2464      !
2465   END SUBROUTINE s_sh94
2466
2467
2468   SUBROUTINE s_sf12
2469      !!----------------------------------------------------------------------
2470      !!                  ***  ROUTINE s_sf12 ***
2471      !!                     
2472      !! ** Purpose :   stretch the s-coordinate system
2473      !!
2474      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
2475      !!                mixed S/sigma/Z coordinate
2476      !!
2477      !!                This method allows the maintenance of fixed surface and or
2478      !!                bottom cell resolutions (cf. geopotential coordinates)
2479      !!                within an analytically derived stretched S-coordinate framework.
2480      !!
2481      !!
2482      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
2483      !!----------------------------------------------------------------------
2484      INTEGER  ::   ji, jj, jk           ! dummy loop argument
2485      REAL(wp) ::   zsmth               ! smoothing around critical depth
2486      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
2487      REAL(wp) ::   ztmpu, ztmpv, ztmpf
2488      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1
2489      !
2490      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
2491      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
2492      !!----------------------------------------------------------------------
2493      !
2494      CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2495      CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2496
2497      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
2498      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
2499      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
2500      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
2501
2502      DO ji = 1, jpi
2503         DO jj = 1, jpj
2504
2505          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
2506             
2507              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
2508                                                     ! could be changed by users but care must be taken to do so carefully
2509              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
2510           
2511              zzs = rn_zs / hbatt(ji,jj) 
2512             
2513              IF (rn_efold /= 0.0_wp) THEN
2514                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
2515              ELSE
2516                zsmth = 1.0_wp 
2517              ENDIF
2518               
2519              DO jk = 1, jpk
2520                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
2521                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
2522              ENDDO
2523              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
2524              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
2525 
2526          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
2527
2528            DO jk = 1, jpk
2529              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
2530              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
2531            END DO
2532
2533          ELSE  ! shallow water, z coordinates
2534
2535            DO jk = 1, jpk
2536              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2537              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
2538            END DO
2539
2540          ENDIF
2541
2542          DO jk = 1, jpkm1
2543             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
2544             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
2545          END DO
2546          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
2547          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
2548
2549          ! Coefficients for vertical depth as the sum of e3w scale factors
2550          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
2551          DO jk = 2, jpk
2552             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
2553          END DO
2554
2555          DO jk = 1, jpk
2556             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
2557             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
2558             gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
2559          END DO
2560
2561        ENDDO   ! for all jj's
2562      ENDDO    ! for all ji's
2563
2564      DO ji=1,jpi-1
2565        DO jj=1,jpj-1
2566
2567           ! extend to suit for Wetting/Drying case
2568           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj)
2569           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1)
2570           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
2571           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
2572           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
2573           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
2574                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
2575           DO jk = 1, jpk
2576              IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN
2577                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
2578                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
2579              ELSE
2580                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
2581                       &              / ztmpu
2582                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
2583                       &              / ztmpu
2584              END IF
2585
2586              IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN
2587                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
2588                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
2589              ELSE
2590                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
2591                       &              / ztmpv
2592                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
2593                       &              / ztmpv
2594              END IF
2595
2596              IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN
2597                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)                 &
2598                       &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
2599              ELSE
2600                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               &
2601                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               &
2602                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               &
2603                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
2604              END IF
2605
2606             ! Code prior to wetting and drying option (for reference)
2607             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
2608             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
2609             !
2610             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
2611             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) )
2612             !
2613             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
2614             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
2615             !
2616             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
2617             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) )
2618             !
2619             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 &
2620             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 &
2621             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 &
2622             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               &
2623             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
2624
2625             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
2626             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
2627             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
2628             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
2629             !
2630             e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
2631             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
2632             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
2633          END DO
2634
2635        ENDDO
2636      ENDDO
2637      !
2638      CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
2639      CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
2640      CALL lbc_lnk(e3w_0 ,'T',1.)
2641      CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
2642      !
2643      CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      )
2644      CALL wrk_dealloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
2645      !
2646   END SUBROUTINE s_sf12
2647
2648
2649   SUBROUTINE s_tanh()
2650      !!----------------------------------------------------------------------
2651      !!                  ***  ROUTINE s_tanh***
2652      !!                     
2653      !! ** Purpose :   stretch the s-coordinate system
2654      !!
2655      !! ** Method  :   s-coordinate stretch
2656      !!
2657      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
2658      !!----------------------------------------------------------------------
2659      INTEGER  ::   ji, jj, jk       ! dummy loop argument
2660      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
2661      REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
2662      REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
2663      !!----------------------------------------------------------------------
2664
2665      CALL wrk_alloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
2666      CALL wrk_alloc( jpk,   z_esigt, z_esigw )
2667
2668      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
2669      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
2670
2671      DO jk = 1, jpk
2672        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
2673        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
2674      END DO
2675      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
2676      !
2677      ! Coefficients for vertical scale factors at w-, t- levels
2678!!gm bug :  define it from analytical function, not like juste bellow....
2679!!gm        or betteroffer the 2 possibilities....
2680      DO jk = 1, jpkm1
2681         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
2682         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
2683      END DO
2684      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
2685      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
2686      !
2687      ! Coefficients for vertical depth as the sum of e3w scale factors
2688      z_gsi3w(1) = 0.5_wp * z_esigw(1)
2689      DO jk = 2, jpk
2690         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
2691      END DO
2692!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
2693      DO jk = 1, jpk
2694         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2695         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
2696         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
2697         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
2698         gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
2699      END DO
2700!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
2701      DO jj = 1, jpj
2702         DO ji = 1, jpi
2703            DO jk = 1, jpk
2704              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2705              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2706              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2707              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2708              !
2709              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2710              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2711              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2712            END DO
2713         END DO
2714      END DO
2715      !
2716      CALL wrk_dealloc( jpk,   z_gsigw, z_gsigt, z_gsi3w )
2717      CALL wrk_dealloc( jpk,   z_esigt, z_esigw          )
2718      !
2719   END SUBROUTINE s_tanh
2720
2721
2722   FUNCTION fssig( pk ) RESULT( pf )
2723      !!----------------------------------------------------------------------
2724      !!                 ***  ROUTINE fssig ***
2725      !!       
2726      !! ** Purpose :   provide the analytical function in s-coordinate
2727      !!         
2728      !! ** Method  :   the function provide the non-dimensional position of
2729      !!                T and W (i.e. between 0 and 1)
2730      !!                T-points at integer values (between 1 and jpk)
2731      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2732      !!----------------------------------------------------------------------
2733      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
2734      REAL(wp)             ::   pf   ! sigma value
2735      !!----------------------------------------------------------------------
2736      !
2737      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
2738         &     - TANH( rn_thetb * rn_theta                                )  )   &
2739         & * (   COSH( rn_theta                           )                      &
2740         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
2741         & / ( 2._wp * SINH( rn_theta ) )
2742      !
2743   END FUNCTION fssig
2744
2745
2746   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2747      !!----------------------------------------------------------------------
2748      !!                 ***  ROUTINE fssig1 ***
2749      !!
2750      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
2751      !!
2752      !! ** Method  :   the function provides the non-dimensional position of
2753      !!                T and W (i.e. between 0 and 1)
2754      !!                T-points at integer values (between 1 and jpk)
2755      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2756      !!----------------------------------------------------------------------
2757      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
2758      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
2759      REAL(wp)             ::   pf1   ! sigma value
2760      !!----------------------------------------------------------------------
2761      !
2762      IF ( rn_theta == 0 ) then      ! uniform sigma
2763         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2764      ELSE                        ! stretched sigma
2765         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
2766            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
2767            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
2768      ENDIF
2769      !
2770   END FUNCTION fssig1
2771
2772
2773   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2774      !!----------------------------------------------------------------------
2775      !!                 ***  ROUTINE fgamma  ***
2776      !!
2777      !! ** Purpose :   provide analytical function for the s-coordinate
2778      !!
2779      !! ** Method  :   the function provides the non-dimensional position of
2780      !!                T and W (i.e. between 0 and 1)
2781      !!                T-points at integer values (between 1 and jpk)
2782      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2783      !!
2784      !!                This method allows the maintenance of fixed surface and or
2785      !!                bottom cell resolutions (cf. geopotential coordinates)
2786      !!                within an analytically derived stretched S-coordinate framework.
2787      !!
2788      !! Reference  :   Siddorn and Furner, in prep
2789      !!----------------------------------------------------------------------
2790      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
2791      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
2792      REAL(wp), INTENT(in   ) ::   pzb            ! Bottom box depth
2793      REAL(wp), INTENT(in   ) ::   pzs            ! surface box depth
2794      REAL(wp), INTENT(in   ) ::   psmth          ! Smoothing parameter
2795      !
2796      INTEGER  ::   jk             ! dummy loop index
2797      REAL(wp) ::   za1,za2,za3    ! local scalar
2798      REAL(wp) ::   zn1,zn2        !   -      -
2799      REAL(wp) ::   za,zb,zx       !   -      -
2800      !!----------------------------------------------------------------------
2801      !
2802      zn1  =  1._wp / REAL( jpkm1, wp )
2803      zn2  =  1._wp -  zn1
2804      !
2805      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
2806      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
2807      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
2808      !
2809      za = pzb - za3*(pzs-za1)-za2
2810      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
2811      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
2812      zx = 1.0_wp-za/2.0_wp-zb
2813      !
2814      DO jk = 1, jpk
2815         p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
2816            &          zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
2817            &               (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
2818        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
2819      END DO
2820      !
2821   END FUNCTION fgamma
2822
2823   !!======================================================================
2824END MODULE domzgr
Note: See TracBrowser for help on using the repository browser.