source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 5621

Last change on this file since 5621 was 5621, checked in by mathiot, 5 years ago

UKMO_ISF: upgrade to NEMO_3.6_STABLE (r5554)

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