source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 4726

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

ISF branch: change name of 2 variables (icedep ⇒ risfdep and lmask ⇒ ssmask), cosmetic changes and add ldfslp key

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